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CHAPTER  I.  INTRODUCTION 


The  work  presented  in  this  thesis  is  concerned  with  modeling  of 
millimeter-wave  impact  avalanche  ;transit-;time  (IMPATT)  diodes.  The  I^TATT 
is  currently  the  most  important  semiconductor  device  for  generation 
and  amplification  of  power  in  the  millimeter-wave  frequency  range. 
Difficulties  arise  in  the  modeling  of  these  devices  because  of  their 
submicron  dimensions  and  high  operating  frequencies.  Device  models 
based  on  the  well  known  "drift  and  diffusion"  description  of 
electron  and  hole  transport  are  not  expected  to  always  be  applicable 
for  such  small  dimensions  and  high  frequencies,  so  it  is  necessary  to 
have  more  detailed  models,  both  to  establish  the  circumstances  under 
which  the  drift  and  diffusion  approximation  fails  and  to  model  devices 
to  which  drift  and  diffusion  based  analyses  cannot  be  applied.  This 
thesis  describes  the  development  of  a  more  general  transport  model, 
the  development  of  a  computer  simulation  based  on  the  model,  and  the 
results  obtained  using  the  simulation. 

This  introductory  chapter  consists  of  three  sections.  The  first 
two  establish  the  context  of  the  present  work  by  providing  reviews  of 
the  principles  and  performance  capabilities  of  IMPATT  diodes  and  of 
approaches  to  modeling  charge  transport  in  semiconductors.  The  third 
section  describes  the  organization  of  the  remainder  of  the  thesis. 

1.1  IMPATT  Devices 

1.1.1  Basic  Operating  Principles  and  Experimental  State  of 
the  Art.  The  use  of  carrier  transit-time  effects  to  produce  negative 
resistance  in  a  semiconductor  device  was  first  proposed  by  Shockley 


in  195^.^  In  1953,  Read^  presented  an  analysis  of  a  diode  structure 
in  which  use  of  transit-time  effects  would  he  combined  with  carrier 
injection  by  impact  ionization.  The  fundamental  characteristics  of 
Read-type  device  operation  can  be  understood  by  considering  the 
operation  of  an  idealized  device,  whose  doping  and  dc  electric  field 
profiles  are  shown  in  Fig.  1.1.^  The  figure  also  shows  a  sinusoidal 
RF  terminal  voltage  and  the  idealized  forms  of  the  injected  and  induced 
current  wave  forms  which  result  when  the  magnitude  of  the  dc  bias 
voltage  is  just  below  that  required  for  reverse  breakdown. 

The  doping  profile  is  chosen  so  that  under  reverse  bias  a 
narrow  "ionization"  region  of  high  field  exists  near  the  p-n  junction. 
At  the  beginning  of  the  RF  cycle,  the  field  in  this  region  is  slightly 
below  threshold  for  avalanche  breakdown  by  impact  ionization.  As 
terminal  voltage  increases,  the  threshold  is  passed,  and  the  number  of 
carriers  in  the  ionization  region  begins  to  increase.  This  continues 
until  the  midpoint  of  the  cycle  when  the  field  in  the  ionization  region 
drops  once  sigain  below  threshold.  The  electrons  generated  in  the 
ionization  region  are  injected  into  the  lower-field  "drift"  region, 
where,  if  the  field  is  strong  enough,  they  travel  at  an  approximately 
constant  velocity.  If  the  length  of  the  drift  region  is  such  that  the 
drift  transit  time  is  half  the  RF  period,  the  motion  of  the  electrons 
Induces  a  flow  of  current  in  the  diode  terminals  which  is  approximately 
l80  degrees  out  of  phase  with  the  terminal  voltage,  giving  rise  to 
negative  resistance.  As  the  cycle  ends,  the  drifting  electrons  are 
collected  by  the  rl^t-hand  contact,  and  the  process  repeats.  All 
devices  which  operate  by  this  combination  of  injection  by  Impact 
ionization  and  drift  across  a  depleted  region  are  known  as  IMPATT 
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diodes,  though  their  structure  may  differ  considerably  from  that 
shown  in  Fig.  1.1. 

In  the  operation  of  a  real  IMPATT  device,  various  departiires 
from  the  ideeJ.  behavior  summarized  in  Fig.  1.1  occur.  The  space 
charge  of  the  carriers  which  accumulate  during  the  first  part  of  the  RF 
cycle  depresses  the  field  in  the  ionization  region  so  that  injection 
occurs  before  the  midpoint  of  the  cycle.  This  reduces  the  lag  between 
terminal  voltage  and  current  and  degrades  the  device  efficiency.  Other 
effects  which  lower  efficiency  include  diffusive  spreading  of  the 
injected  carrier  pulse,  impact  ionization  in  the  drift  region,  and 
carrier  drift  at  nonsatxirated  velocities. 

Because  of  difficulties  with  device  fabrication,  IMPATT  mode 
operation  was  not  realized  until  1965,  when  Lee  et  al.**  succeeded  in 
obtaining  the  first  oscillations  from  a  Read-type  diode.  About  the 
same  time,  Johnson  et  al.®  obtained  oscillations  from  a  simpler  p-n 
diode  structure.  IMPATT  fabrication  and  circuit  technology  have 
developed  steadily  since  then,  and  today  IMPATTs  are  widely  used  as 
sources  and  amplifiers  in  low  and  medium  power  microwave  and 
millimeter-wave  systems.  Silicon  IMPATTs  in  particular  are  presently 
the  most  important  solid  state  power  source  at  high  microwave  and 
millimeter-wave  frequencies.®  The  experimental  power  and  frequency 
state  of  the  art  for  a  variety  of  semiconductor  devices  is  shown  in 
Fig.  1.2.’"^*  This  shows  that  Si  IMPATTs  are  currently  the  highest 
power  semiconductor  devices  for  millimeter-wave  power  generation, 
and  that  their  useful  frequency  range  extends  to  several  hvmdred 
gigahertz. 
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1.1.2  Status  of  IMFATT  Modeling.  IMPATT  operation  is 
inherently  nonlinear,  so  detailed  modeling  of  large-signal  operation 
requires  the  use  of  numerical  methods.  Most  IMPATT  models  to  date  have 
been  intended  for  microwave  devices  and  have  been  based  on  the  conven¬ 
tional  drift  and  diffusion  description  of  carrier  transport.  In  his 
original  paper,  Bead^  assumed  saturated  drift,  with  equal  velocities  and 
ionization  rates  for  electrons  and  holes,  and  ignored  diffusion.  Other 
workers  have  developed  small-signal  models  which  allowed  for  unequal 
saturated  velocities  and  ionization  rates, arbitrary  doping  profiles, 
field-dependent  velocities, and  diffusion. Gilden  and  Hines^** 
derived  a  useful  small-signal  equivalent  circuit  for  Read-type  IMPATTs, 
showing  the  tuning  effects  of  the  dc  bias  current.  Evans  and  Haddad^® 
developed  the  first  closed-form  expression  for  IMPATT  large-signal 
Impedance,  using  the  assumption  of  a  small  phase  angle  associated  with 
the  drift  transit  time.  This  assumption,  together  with  that  of  saturated 
drift,  was  removed  in  a  large-signal  model  developed  by  Greiling  and 
Haddad.^®  A  well  known  finite-difference  simulation  based  on  a  compara¬ 
tively  complete  version  of  the  drift-diffusion  model  was  developed  by 
Scharfetter  and  Gummel.®^  Subsequent  workers  have  greatly  improved 
the  efficiency  of  finite-difference  simulations  based  on  the  drift- 
diffusion  model,  and  have  reached  a  better  understanding  of  the 
numerical  diffusion  which  is  associated  with  finite-difference 
schemes . ® ^  ® 

A  few  attempts  have  been  made  to  develop  IMPATT  models  which 
accoimt  for  additional  physical  effects.*’”"*^  It  will  be  shown  that 
none  of  these  represents  a  self-consistent  treatment  using  a 
transport  model  better  than  conventional  drift  and  diffusion. 
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Development  of  such  a  model  is  one  focus  of  the  present  work,  and  a 
brief  review  of  the  spectrvun  of  possible  transport  models  will  there¬ 
fore  now  be  given  to  establish  a  context  for  the  work  described  in 
subsequent  chapters. 

1.2  Models  of  Electron  Transport  in  Semiconductors 

1.2.1  A  Hierarchy  of  Approaches.  The  range  of  treatments  of 
electron  transport  in  semiconductors  is  conveniently  described  with 
reference  to  the  hierarchy  of  approaches  shown  in  Fig.  l.B.**^  Electron 
transport  in  semiconductors  is  fundamentally  quantum  mechanical  in 
nature  (because  the  deBroglie  wavelength  of  electrons  is  not  small 
with  respect  to  interatomic  spacings),  but  can  often  be  modeled  using 
the  quasi-free-particle  (QFP)  approximation  shown  in  the  center  of  the 
figure.  In  the  QFP  approximation,  individual  electrons  are  treated  as 
classical  particles  with  effective  mass  supplied  by  band  theory.  The 
concept  of  the  hole  is  used  to  describe  charge  transport  due  to  empty 
states  in  the  valence  band.  The  approximation  is  based  on  the  assumption  that 
collisions  and  acceleration  due  to  externally  applied  electric  fields 
can  be  treated  as  perturbations  to  the  band  structure  description  of 
the  perfect  lattice.  This  is  generally  true  if  the  field  is  not  so 
large  as  to  invalidate  the  use  of  Bloch  functions  for  the  electron 
states,  and  if  characteristic  distances  other  than  interatomic 
spacings  are  large  compared  to  the  size  of  an  electron  wave  packet."*^ 
Once  the  band  structure  and  scattering  rates  are  established,  QFP 
modeling  of  carrier  transport  can  be  accomplished  by  use  of  Monte 
Carlo  or  Rees  iterative'*'*  techniques  to  solve  the  phase-space 
transport  equation'**  with  the  appropriate  collision  term.  However, 
such  methods  are  generally  too  expensive  for  routine  application 
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1.0  l.ime  domain  simulation  of  devices,  altliough  they  are  used  extensively 
’’or  modeling  carrier  transport  phenomena. 

When  certain  simplifications  are  made  in  order  to  produce  more 
economical  device  simulations,  sub-QFP  models  result.  Intead  of 
working  in  terms  of  individual  carriers  or  the  exact  form  of  the  dis¬ 
tribution  function,  these  models  use  approximate  descriptions  of  the 
distribution  function.  The  distribution  may  be  assumed  to  be  in  some 
parameterized  form,  such  as  displaced  Maxwellian  or  a  truncated  series 
of  Legendre  polynomials,  after  which  equations  for  the  unknown  param¬ 
eters  can  be  found  from  the  phase-space  transport  equation.  Alterna¬ 
tively,  Monte  Carlo  data  can  be  used  to  evaluate  the  phase-space 
collision  term  directly,  and  the  results  used  in  simple  energy  and 
momentum  balance  relations.  This  latter  approach  requires  no  assump¬ 
tions  about  the  form  of  the  velocity  distribution,  but  intuitively 
chosen  balance  relations  are  in  fact  inconsistent  with  the  phase- 
space  equation  in  spatially  inhomogeneous  situations. 

The  transport  picture  can  be  further  simplified  by  assuming  that 
the  velocity  distribution  is  always  in  equilibrium  with  the  local  elec¬ 
tric  field.  These  will  be  referred  to  as  "static"  conditions  elsewhere 
in  the  thesis.  The  assumption  is  valid  if  the  change  in  field  strength 
seen  by  moving  carriers  is  small  during  the  time  required  to  reach 
equilibrium.  In  this  "static"  approximation,  carrier  motion  is 
commonly  described  using  the  conventional  drift-diffusion  equation  with 
field-dependent  drift  velocity  (or  mobility),  diffusion  coefficients, 
and  ionization  rates.  In  some  situations  where  the  static  model  is 
almost  appropriate  for  describing  carrier  transport,  additional  effects 
can  be  accommodated  by  adding  extra  transport  parameters.  A 
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"reparameterized"  model  has  been  used,  for  example,  in  dead  space 

^  6 

modeling  of  the  impact  ionization  process. 

Above  the  QFP  approximation  in  Fig.  1.3  are  the  partial  and  full 
quantum  regimes.  In  the  partial  quantum  regime,  the  QFP  approximation 
is  substantially  valid,  but  additional  quantum  effects,  such  as  band-to- 
band  tunneling  or  size  quantization  in  one  spatial  dimension,  are  im¬ 
posed  on  the  basic  framework.  When  characteristic  dimension  becomes 
extremely  small,  full  quantum  treatment  of  transport  becomes  necessary. 
This  is  still  primarily  the  realm  of  the  theoretical  physicist. 

1.2.2  Applications  to  IMPATT  Modeling.  The  preceeding  discussion 
illuminates  why  the  applicability  of  the  static  drift-diffusion  model 
to  IMPATT  modeling  becomes  suspect  in  the  case  of  devices  operating  at 
millimeter-wave  frequencies.  The  drift-diffusion  mod.-'f  .'.ssumct  vjijuili- 
brium  between  the  carrier  velocity  distributions  and  local  electric 
field,  but  in  millimeter-wave  IMPATTa  such  equilibrium  frequently  may 
not  exist  because  it  is  possible  for  carriers  in  these  devices  to 
experience  significant  changes  in  field  strength  during  the  time  re¬ 
quired  to  reach  equilibrium.  There  are  two  main  peasons  why  depar¬ 
tures  from  equilibrium  will  be  more  significant  at  millimeter-wave 
frequencies  than  at  microwave  frequencies.  First,  design  length 
decreases,  and  doping  levels  increase  with  frequency,  so  that  drifting 
carriers  travel  through  steeper  gradients  of  the  field  strength  in 
millimeter-wave  devices  than  in  microwave  devices.  This  will  affect 
even  the  dc  behavior.  Second,  the  maximum  rate  of  change  of  terminal 
voltage,  hence  of  internal  field  strength,  tends  to  increase  with  fre¬ 
quency,  though  this  tendency  is  offset  somewhat  by  the  decrease  in 
RF  amplitude  which  comes  with  decreasing  device  length. 
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Use  of  a  full  QFP  model  vould  give  full  knowledge  of  the 
carrier  velocity  distributions  at  all  times.  As  was  noted,  however, 
such  models  are  generally  expensive  when  applied  to  simulation  of 
devices.  The  reason  for  this  is  that  they  keep  track  of  excessive 
information:  the  position  and  velocity  of  each  individual  carrier.  Much 
of  this  information  is  not  of  first-order  importance  to  a  device  simula¬ 
tion,  since,  for  the  purpose  of  determining  device  terminal  currents 
and  voltages,  only  the  spatial  distributions  of  carrier  concentration 
and  average  velocity  are  re'-iuired. 

One  attractive  approach  to  filling  the  need  for  a  nonstatic  model 
for  millimeter-wave  IMPATTs  is  the  method  of  conservation  of  energy  and 
momentum,  which  falls  in  the  sub-QFP  regime  of  the  hierarchy  shown  in 
Fig.  1.3.  This  method  is  not  based  on  the  static  assumption  of 
carrier  field  equilibrium  and  is  much  more  economical  for  device 
simulation  than  full  QFP  methods.  A  presentation  of  the  equations  of 
energy  and  momentum  conservation  in  transport  in  semiconductors,  to¬ 
gether  with  an  extensive  discussion  of  their  relationship  to  various 
lower  order  models,  has  been  given  by  Blotekjaer. 

The  nature  of  the  energy  and  momentum  conserving  model  for 
carrier  transport  can  be  briefly  described  as  follows.  The  model 
consists  of  Poisson's  equation  for  the  electric  field  gradient, 
together  with  transport  equations  for  three  quantities  as  functions 
of  space  and  time:  carrier  concentration,  average  momentvun,  and 
average  energy.  The  first  two  quantities  are  required  for  determin¬ 
ing  device  terminal  behavior.  The  momentum  equation  keeps  track  of 
the  various  contributions  to  the  momentum  of  the  aggregate  of  carriers, 
such  as  gains  due  to  acceleration  by  the  field  and  losses  due  to 
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collisions.  This  means  that  the  average  carrier  velocity  is  not  a 
static  function  of  field.  The  third  quantity,  average  energy,  is 
not  required  for  calculating  device  voltages  and  currents,  but  gives 
the  width  of  the  velocity  distribution,  which  affects  the  rate  of 
carrier  diffusion  and,  to  a  first  approximation,  determines  the  rates 
of  collisions  which  affect  momentum  and  determines  the  rate  of  impact 
ioni zation. 

Concentration,  average  velocity,  and  average  energy  are  pro¬ 
portional  to  the  zeroth,  first,  and  second  moments  of  the  velocity 
distribution.  In  the  QFP  approximation,  the  distribution  is  governed 
by  the  phase-space  transport  equation,  so  the  carrier,  momentum,  and 
energy  transport  equations  required  under  the  energy  and  momentum  con¬ 
serving  model  can  be  obtained  by  taking  the  first  three  velocity  moments 
of  the  phase-space  equation.  Using  the  resulting  transport  equations 
in  the  energy  and  momentum  conserving  model  guarantees  that  it  will  be 
consistent  with  the  phase-space  equation  to  second  order  in  the 
velocity  coordinate. 

1.3  Outline  of  the  Present  Study 

The  goals  of  this  study  are  to  apply  the  principles  of  energy 
and  momentum  conservation  to  simulation  of  millimeter-wave  Si  IMPATT 
diodes,  to  produce  a  computer  simulation  embodying  the  principles, 
to  use  the  simiilation  to  examine  device  behavior  under  dc  and  large- 
signal  conditions,  and,  by  comparison  with  results  obtained  using  a 
conventional  simulation,  to  establish  the  limits  of  applicability 
of  drift-diffusion  simulation. 
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The  organization  of  the  thesis  is  as  follows.  Chapter  II 
develops  energy  and  momentum  conserving  transport  equations  appropriate 
for  carriers  in  Si.  Chapter  III  examines  numerical  methods  based  on 
finite-difference  approximation  of  the  transport  equations  derived  in 
Chapter  II.  A  stable,  accurate,  and  efficient  numerical 
procedure  for  their  solution  is  developed.  Chapter  IV  presents  results 
from  computer  simulation  of  millimeter-wave  Si  IMPATT  diodes.  Results 
from  the  new  simulation  are  compared  with  results  obtained  using  a 
conventional  drift-diffusion  simulation  for  a  variety  of  device  lengths 
and  operating  frequencies,  and  reasons  for  the  observed  differences 
are  discussed.  The  effects  of  various  energy  boundary  conditions  and 
of  realistic  contact  regions  are  also  examined.  Chapter  V  contains 
discussion,  conclusions,  and  suggestions  for  further  research. 
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CHAPTER  II.  THE  TRANSPORT  MODEL 


This  chapter  describes  the  energy  and  momentum  conserving 
transport  model  which  has  been  developed  to  provide  a  better  des¬ 
cription  of  carrier  transport  in  Si  IMPATT  diodes  than  that  provided 
by  the  conventional  drift-diffusion  model-  The  contents  of  the  chapter 
are  as  follows.  Section  1  presents  a  derivation  of  the  collisionless 
forms  of  the  carrier,  energy,  and  momentum  transport  equations,  and 
discusses  the  physical  interpretations  of  their  various  terms.  In 
Section  2,  functional  forms  and  nvimerical  values  are  obtained  for  terms 
describing  the  effects  of  collision  processes.  Section  3  discusses 
the  differences  between  the  resulting  model  and  other  nonstatic  IMPATT 
models,  and  also  shows  how,  under  certain  conditions,  the  model  limits 
to  the  conventional  one.  Section  U  is  a  general  discussion  and 
summary  of  the  chapter. 

2.1  The  Collisionless  Transport  Equations 

This  section  develops  collisionless  transport  equations  for  carriers, 
carrier  energy  and  carrier  momentum  in  Si.*  The  development  presented 
here  and  throughout  much  of  the  remainder  of  this  chapter  is  given 
in  terms  of  electron  transport  only.  The  extensions  to  hole  transport 
are  readily  apparent. 

2.1.1  Distribution-Independent  Analysis.  The  starting  point 
for  building  the  energy  and  momentum  conserving  transport  model  is 
the  quasi-free-particle  approximation  described  in  Chapter  I. 

*Much  of  the  material  In  this  section  is  based  on  the  treatment  of 
Duderstadt  and  Martin."*® 
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The  approximation  allows  carriers  to  be  treated  as  classical  particles 
with  effective  masses  and  collision  rates  determined  from  the  -nergy 


band  structure  and  perturbation  analysis.  The  motion  of  the  carriers 
can  be  followed  by  keeping  track  of  their  positions  and  velocities.  It 
is  convenient  to  do  so  by  keeping  track  of  the  carrier  distribution 
function,  which  gives  the  carrier  concentration  in  six-dimensional 
(location  and  velocity)  phase  space  as  a  function  of  time. 

In  the  absence  of  collisions,  the  concentration  N(r,v,t)  at  the 
phase  point  (r,v)  and  time  t  will,  according  to  the  Liouville 
theorem,"*®  follow  its  trajectory  in  phase  space  unchanged.  A  short 
time  At  later,  it  will  reach  the  point  [r  +  vAt ,v  +  (F/m)At], 
where  F  is  the  force  acting  and  ra  the  effective  mass  in  the  neighbor¬ 
hood  of  (r,v).  The  only  change  in  N  between  the  two  points  will  be  due 

to  collisions  which  scatter  carriers  into  or  away  from  the  neighborhood 
of  (r,v).  This  can  be  described  by  writing 

Nfr  +  vAt,v  +  (F/m)At.t  +  At]  -  N(r.v,t)  _  "^c^  ,  . 

At  "  At  ’  ^  ^  ^ 


where  6  N  is  the  change  due  to  collisions.  In  the  limit  as  At  approaches 
c 

zero,  Eq.  2.1  becomes 


at 


-  v»VN 


+ 


(2.2) 


'  '  c 

Equation  2.2  is  the  phase-space  transport  equation  which  des¬ 
cribes  the  motion  of  carriers  in  the  QFP  approximation.  The  well 
known  Boltzmann  transport  equation  is  similar  to  Eq.  2.2,  with  a 
collision  term  in  the  particular  form  which  describes  collision 
effects  in  a  dilute  gas. 
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Since  N  is  a  phase-space  density,  the  densities  of  carriers, 
mean  carrier  momentum,  and  mean  carrier  energy  in  coordinate  space 


can  be  defined  in  terms  of  the  zeroth,  first,  and  second  velocity 
moments  of  N: 


and 


n(r,t)  = 
m*n(r,t)u{r,t) 

n(r,t)w(r,t)  = 


I  N(r,v,t)  d*v  , 

=  I  mvN(r,v,t)  d*v 

I  -I  mjvPN(r,v,t)  dV  , 


(2.3) 


(2.5) 


where  n  is  the  carrier,  concentration,  u  is  the  average  velocity,  w  is 
the  average  energy,  and  m*  is  the  average  over  effective  mass.  It 
will  be  assumed  that  m  is  constant.  This  is  reasonable  so  long  as 
most  of  the  carriers  are  fairly  close  to  energy  minima. 

The  energy  and  momentum  conserving  transport  model  consists  of 
transport  equations  for  the  quantities  n,  u,  and  w.  The  equations  can 
be  derived  by  applying  the  method  of  moments  to  Eq.  2.2,  the  phase- 
space  transport  equation,  as  will  now  be  shown.  It  is  convenient  to 
ignore  the  collision  terra  in  Eq.  2.2  for  the  time  being;  its  effect  on 
the  equations  to  be  derived  here  is  taken  up  in  Section  2.2. 

Part  of  the  operation  of  taking  a  velocity  moment  of  Eq.  2.2  can 
be  performed  without  regard  for  the  form  of  a  particular  velocity 
moment  operator.  If  ^  is  any  such  operator,  a  general  moment  of  Eq. 
2.2  is  given  by 


-  vVN 


d’v 


(2.6) 


Since  t|(  is  a  function  of  v  alone,  Eq.  2.6  can  be  written  as 
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Nip  d^v 


=  “  ^  *1  ^  m  */ 


iP  dV  ,  (2.7) 


where  it  is  assumed  that  N  approaches  zero  rapidly  enough  at  the  limits 
of  integration  that  the  quantity  /  7^(Ni|))  d^v  is  always  negligible. 

The  carrier  concentration  is  the  zeroth  velocity  moment  of  N.  With 
=  1,  Eq.  2.7  becomes 


on  n  /  —\ 


(2.8) 


This  is  the  usual  carrier  continuity  equation,  though  no  diffusion  term 
appears  in  it  explicitly.  The  velocity  u  is  the  true  average  over  the 
velocity  distribution  rather  than  the  field-dependent  "drift  velocity" 
used  in  the  static  drift-diffusion  model,  so  that  the  diffusion  effects 
which  must  be  treated  there  in  a  separate  term  are  here  incorporated 


into  u. 


For  velocity,  tj/  equals  v,  and  Eq.  2.7  becomes 


-  V*  vvN  d^v  + 


-•f  N[y  X.  d’v 

ml  V  1  3v.  I 
;  '•1 

=  -  V »!  w  N  d^v  +  N[i1  d^v  ,  (2.9) 


where  [ll  is  the  identity  tensor.  The  left-hand  side  of  Eq.  2.9  can  be 
expanded  and  rewritten  using  Eq.  2.8,  resulting  in 

^  =  -  —  V»  wN  d*v  +  —  +  —  V*(nu)  .  (2.10) 

3x  n  in  n 

The  remaining  integral  can  be  rewritten  as 

V«  wN  d^v  =  n^*^|  “  u)(v  -  u)K  d^v  +  |  (uv  +  vu  -  uu)N  d^-i^ 

=  ^  V-[i  [P]  +  min  ,  (2.11) 
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where  [P]  is  a  tensor  defined  by 


P.  ,  =  m 

ij 


I  N(v^  -  -  Uj)  d^v  . 


(2.12) 

The  remainder  of  the  right-hand  side  of  Eq.  2,11  can  be  expanded  to 


—  7*(uun) 
n 


1  V  4.  ilL-l  =  i 

n  f  ax.  ax.)  n 


=  —  u7*{nu)  +  u*Vu 

n 

Substituting  for  the  integral  in  Eq.  2.10  results  in 


(2.13) 


=  _  u*Vu  +  -  -  —  7.[P] 

at  m  nm 


(2.11+) 


Equation  2.ll*  describes  the  transport  of  average  carrier  velocity, 
in  which  the  divergence  of  [P]  plays  the  role  of  a  force  field  which 
contributes  to  acceleration.  The  first  term  on  the  right  of  Eq. 

2.1I+  is  related  to  the  divergence  of  velocity  flux.  It  is  not  in 
"conservation  form,"**®  where  the  derivative  operator  acts  on  the 
entire  term,  as  in  the  corresponding  term  in  Eq.  2.8.  This  is  because 
u  is  not,  in  fact,  a  conserved  quantity  in  carrier  transport;  it 
represents  average,  rather  than  total,  momentum,  and  it  is  total 
moment+ira  which  is  conserved  physically. 

For  energy,  1)1  equals  mvv/2.  Substitution  in  Eq.  2.7  gives 

^  7*1  vNv«v  d®v  +  ■— *1  N7^(v*v)  d®v  .  (2.15) 

The  second  of  the  two  integral  terms  in  Eq.  2.15  can  be  simplified; 

1*1  N7^(vv)  d®v  =  ^  I  (vp  d®v 

=  F*|  Nv  d®v 

=  nF*u  .  (2.16) 
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( 

* 


Since  /  Nu*u(v  -  u)  d^v  is  zero,  the  first  of  the  integral  terms 


can  he  written  as 


^  V*  I  vNvv  d^v  =  ^  7.  |(vv  -  u*u)(v  -  u)N  d^v  +  | 


uv«vN  d*v 


^7*11  [ (v  -  u)»(v  -  u)  +  2u*(v  -  u)]Tv-u)N  d^ 


^2  - 
+  —  wnu 
m 


=  7*(u*[P]  +  q)  +  nu»Vw  +  wV*(nu)  ,  (2.17) 


where  the  vector  q  is  defined  by 


‘^i  “  j  2  ^^^i  “  -  u)*(v  -  u)  d*v  . 


(2.18) 


With  the  use  of  Eqs.  2.8,  2.16  and  2.17,  Eq.  2.15  becomes 


=  -  u*Vw  +  F.u  -  -  v(u*lP]  +  q)  .  (2.19) 

dt  n 


Equation  2.19  describes  the  transport  of  average  energy  w. 
Here  again,  as  in  Eq.  2.lU,  the  flux  divergence  term  u»7w  is 
in  nonconservative  form  because  average  energy  is  not  conserved 
physically.  The  vector  q  is  known  in  fluid  mechanics  as  the  heat 
flow  vector. 
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2.1.2  Dependence  on  the  Velocity  Distribution .  Equations 
2.8,  2.1h^  and  2.19  describe  the  transport  of  velocity-averaged  quan¬ 
tities,  but  their  solution  requires  knowledge  of  the  form  of  the 
velocity  distribution  N,  which  appears  in  the  definitions  of  [P]  and 
q.  Some  form  of  the  velocity  distribution  must  be  assumed  for  the 
purpose  of  evaluating  [P]  and  q.  It  should  contain  n,  u,  and  w  as 
parameters,  so  that  [P]  and  q  will  be  consistent  with  these.  At  the 
same  time,  specifying  the  distribution  beyond  its  second  velocity 
moment  is  impractical,  since  the  energy/  and  momentum  conserving  model 
gives  no  information  as  to  hov  the  higher  moments  change  with  time. 
For  the  purpose  of  finding  [P]  and  q,  it  will  be  assumed  that  the 
velocity  distribution  is  displaced  Maxwellian  in  form,  which 
makes  the  integrations  in  Eqs.  2.12  and  2.l8  particularly  simple. 

The  displaced  Maxwellian  form  whose  first  three  moments  are  consis¬ 
tent  with  n,  u,  and  w  is 


f  m 

m!v  -  u 1 

r  J 

(2.20) 


where  Boltzmann's  constant,  and  is  the  carrier  temperature, 

defined  by 


I  V=  =  H  5  I’'  ■  "I'" 

=  I  ^  (vv  +  u*u  -  2u*v)N  d®v 

=  w  -  ^  u*u  .  (2.21) 

Thus,  the  thermal  component  of  w  is  3k^T^/2,  in  keeping  with  the  usual 
idea  of  temperature. 
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Changing  variables  to  v'  =  v  -  u  simplifies  the  integrations 
in  Eqs.  2.12  and  2.18.  Since  N  is  even  in  v' ,  integration  over  all 
of  a  product  of  N  with  a  function  that  is  odd  in  any  component  of 
v'  will  give  zero.  Equations  2.12  and  2.l8  then  become 


P.  .  =  Nviv'm  d^v' 

ij  j  1  >5 

=  0  ,  i  /  j 

=  ,  i  =  J  (2.22) 

and 


=  0  ,  all  i  .  (2.23) 


Then  Eqs.  2.lU  and  2.19  become 


n  =  - ^  ■  I  (2.2M 

and 

(2.25) 

where  q  is  the  electronic  charge  and  E  is  the  electric  field. 

The  forms  in  Eqs.  2.8,  2.2U,  and  2.25  of  the  collisionless 
transport  equations  are  equivalen’  to  the  forms  of  the  equations  of 
hydrodynamics  describing  the  motion  of  ein  inviscid,  compressible 
fluid.®®  The  terms  arising  from  [P]  will  be  referred  to  as  pressure 
terms  because  the  diagonal  entries  in  [P]  are  Just  the  pressure  of 
an  ideal  gas  with  concentration  n  at  temperature  T^.  The  pressure 
terms  account  for  the  momentum  and  energy  imparted  to  carriers 
traveling  down  a  pressure  gradient  in  the  carrier  "gas." 


3w 

at 


=  -  u»Vw  + 


—  2  r-.  1  — ; 

qE-u  -  V  nu(w  -  —  mu»uj 
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The  energy  and  momentum  conserving  transport  model  is  not 
exact  in  the  QFP  approximation  because  it  is  in  terms  of  only  the 
first  three  velocity  moments  of  N,  rather  than  the  exact  form  of  N. 

It  should  be  noted  that  q  has  been  eliminated  going  from  Eqs.  2.19 
to  2. 2U  because  the  assumed  displaced  Maxwellian  distribution  is 
symmetrical  about  u.  BloteRjaer"’  has  suggested  that  the  model  can 
be  improved  by  retaining  nonzero  q  in  the  energy  transport  equation. 
However,  Bosch  and  Thim*^  used  an  energy  and  momentum  conserving  model 
to  simulate  transferred  electron  devices  and  found  that  the  inclusion 
of  nonzero  q  had  little  effect  on  predicted  device  behavior.  It  will 
be  assumed  here  that  the  transport  model  will  be  sufficiently  accurate 
for  use  in  device  modeling  without  the  inclusion  of  q  in  Eq.  2.25* 

It  will  also  be  assumed  that  carrier  transport  in  IMPATTS  is  one¬ 
dimensional,  so  that  henceforward  the  transport  equations  can  be 
written  in  scalar  form. 

2.2  Collision  Terms 

This  section  develops  terms  for  inclusion  in  the  carrier,  energy, 
and  mcxnentum  transport  equations  which  take  into  account  the  effects  of 
the  phase-space  collision  term.  The  procedure  followed  in  the  previous 
section  was  to  find  velocity  moments  of  the  collisionless  phase- 
space  equation.  The  collision  terms  will  be  found  using  a  different 
procedure,  for  reasons  which  will  be  discussed. 

2.2.1  Physical  Considerations .  The  modifications  which  are  to 
be  made  to  the  carrier,  momentum,  and  energy  tremsport  equations, 

Eqs.  2.8,  2.2k,  and  2.25,  must  account  for  the  rates  at  which  carrier 
concentration,  mean  velocity  and  mean  energy  are  changing  because  of 
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collision  processes.  These  rates  of  change  depend  in  reality  on  the 
exact  form  of  the  distribution  function,  but  this  is  not  knovm  under 
the  energy  and  momentum  conserving  model.  The  problem  becomes  one  of 
accurately  approximating  the  rates  of  change  with  functions  of 
n,  u,  and  w. 

An  obvious  way  of  doing  so  is  to  find  the  first  three  velocity 
moments  of  the  phase-space  collision  term  and  include  them  at  the 
appropriate  points  in  the  analysis  presented  in  the  previous  section. 
This  is  difficult  to  do,  however,  because  the  phase-space  term  can 
seldom  be  known  with  much  accuracy.  It  contains  the  distribution 
function  and  functions  describing  the  average  rates  and  effects  of 
the  various  collision  types;  e.g.,  phonon  emission  and  absorption, 
impact  ionization,  impurity  and  defect  scattering,  and  none  of  these 
is  always  known  exactly.  The  rate  function  for  impact  ionization  is 
.ery  difficult  to  determine,®^"®'*  and  while  functions  are  known  for 
a  number  of  other  collision  processes,®®  these  contain  adjustable 
parameters,  whose  values  are  in  practice  chosen  to  make  theoretical 
pi cdictions  agree  with  experimental  measurements. 

In  effect,  measured  data  are  what  determine  numerical  values 
of  the  phase-space  collision  term.  This  suggests  a  more  direct  way 
of  arriving  at  the  carrier,  energy,  and  momentum  collision  terms: 
that  they  be  evaluated  directly  from  this  same  data.  This  is  the 
procedure  that  will  be  followed  here. 

2.2.2  Forms  of  the  Collision  Terms .  Before  the  collision  terms 
can  be  evaluated,  the  forms  of  their  dependencies  on  the  variables  n, 
u,  and  w  of  the  transport  model  must  be  established.  The  collision 
types  to  be  accounted  for  (considered  here  in  terms  of  their  effects 
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on  electrons )  can  be  divided  into  three  categories  according  to  the 
types  of  carriers  (electrons  and/or  holes)  they  involve: 

1.  Collisions  undergone  by  electrons  which  change  the 
concentration,  energy,  or  momenttun  of  electrons. 

2.  Collisions  undergone  by  holes  which  change  the  concen¬ 
tration,  energy,  or  momentum  of  electrons. 

3.  Electron-hole  interactions  which  change  the  concentra¬ 
tion,  energy  or  momentum  of  electrons. 

Electron-electron  interactions  are  not  considered  because  they  have  no 
effect  on  electron  concentration,  average  energy,  or  average  momentim. 
The  effect  of  electron-electron  scattering  is  in  any  case  to  make  the 
distribution  more  nearly  a  displaced  Maxwellian.®® 

It  will  be  assumed  that  all  types  of  collisions  which  fit  in 
categories  (1)  or  (2)  above  involve  single  carriers,  and  that  their 
per  carrier  rates  are  functions  only  of  average  carrier  energy. 
Blotekjaer  and  Lunde®*  have  shown  that  the  latter  assiamption  is  reason¬ 
able  in  the  case  of  a  displaced  !iEixwellian  distribution.  They  derived 
formulae  for  the  energy  and  momentum  relaxation  times  (defined  below) 
by  taking  moments  of  the  phase-space  collision  term.  Their  results, 
when  written  in  terms  of  w,  are  independent  of  u  to  second  order  in  u. 

2. 2. 2.1  Same-Carrier  Collisions.  Category  (1)  above 
Includes  all  lattice  collisions  undergone  by  electrons.  Of  these,  the 
carrier  concentration  is  affected  by  impact  ionization,  which  contri¬ 
butes  a  generation  rate  to  the  carrier  transport  equation: 

(2.26) 

where  a  is  the  average  per  electron  ionization  rate  euid  is  taken  to  be 


is. 

3t 


=  an 
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a  function  of  w.  Its  dimensions  are  inverse  time,  so  it  is  different 
from  the  conventional  electron  ionization  rate  employed  in  the  static 
model.  The  latter  has  dimensions  of  inverse  length,  as  a  carryover 
from  studies  of  gas  discharges.®^  In  fact,  the  probability  that  an 
individual  carrier  will  cause  ionization  in  unit  time  is  dependent 
upon  its  position  in  the  energy  band  structure,  a  position  which  is, 
in  the  nonstatic  case,  largely  independent  of  its  velocity.  The 
average  of  the  probability  over  a  group  of  carriers  is  therefore, 
in  general,  almost  independent  of  their  average  velocity,  so  the  per- 
unit-time  ionization  rate  a  as  defined  in  Eq.  2.26  is  more  fundamental 
than  the  conventional,  per-unit-di stance  rate. 

Average  momentum  tends  to  be  reduced  by  collision  processes  in 
Category  (l),  since  they  tend  to  randomize  the  velocity  distribution. 
While  the  per-carrier  rates  of  the  processes  are  functions  of  energy, 
the  resulting  rate  of  loss  of  velocity  cannot  be  a  function  of  energy 
alone.  This  can  be  seen  if  two  situations  of  equal  energy  and  concen¬ 
tration  are  considered,  one  in  which  u  is  large  and  one  in  which  u  is 
zero.  If  the  collision  rates  are  functions  only  of  energy,  exactly  the 
same  number  and  types  of  collisions  will  be  taking  place  in  both 
situations.  In  the  first,  u  is  decreasing  comparatively  rapidly  because 
of  the  randomizing  effect  of  collisions.  In  the  second,  u  is  not 
changing  at  all,  since  it  is  already  zero.  The  rate  of  change  of 
velocity  due  to  collisions  is  very  different  from  one  situation  to  the 
other. 

The  velocity  and  energy  dependence  of  momentum  loss  in  collisions 
belonging  to  the  first  category  will  be  accounted  for  by  writing 
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where  is  an  energy-dependent,  effective  Tnomentum  relaxation  time. 

While  the  rate  of  loss  of  average  energy  to  these  collisions  is 
a  function  of  energy  alone,  it  is  convenient  for  the  purposes  of  the 
finite-difference  approximations  to  the  transport  equations  developed 
in  Chapter  III  to  model  the  energy  loss  rate  in  terms  Of  an  energy 
relaxation  time.  The  minimum  energy,  instead  of  being  zero,  is  the 
thermal  energy  Wq  associated  with  the  temperature  of  the  lattice,  so  the 


rate  of  energy  loss  will  be  written  as 
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(2.28) 


where  t  is  an  energy  dependent,  effective  energy  relaxation  time.  The 
w 

use  of  effective  relaxation  times  in  Eqs.  2.27  and  2.28  is  in  accordance 
with  the  forms  of  the  energy  and  momentum  transport  equations  given  by 
Blotekjaer. 


2, 2. 2. 2  Opposite- Carrier  Collisions.  The  second  collision 
category  consists  (from  the  point  of  view  of  electrons)  of  impact 
ionization  by  holes.  These  contribute  to  the  carrier  generation  rate. 


=  3P  ,  (2.29) 

where  8  is  the  per-\uiit-time  hole  ionization  rate  and  depends  on  hole 
energy.  Hole  ionizations  also  affect  electron  average  velocity  and 
energy.  If  it  is  assumed  that  the  electrons  created  by  hole  ioniza¬ 
tions  have  random  velocities,  their  contribution  to  the  total  electron 
momentum  is  zero: 
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Then,  combining  Eqs.  2.29  and  2.30  resxilts  in 


2.2.2. 3  Electron-Hole  Interactions .  The  third  collision 
category  includes  Auger  recombination  and  direct  exchange  of  energy 
and  momentum  between  electrons  and  ho;..es.  Sze*®  gives  an  estimate  for 
the  Auger  recombination  lifetime  which  is  large  compared  to  a  millimeter- 
wave  period  of  oscillation  under  any  circumstances  which  might  normally 
occur  in  an  operating  IMPATT.  Blotekjaer  and  Lunde®®  have  estimated 
the  rates  of  energy  and  momentum  transfer  between  electrons  and  holes 
under  the  assumption  of  a  displaced  liaxwellian  distribution.  Under 
ordinary  conditions,  their  estimates  are  small  in  comparison  to  the 
rates  of  energy  and  momentum  loss  which  result  from  the  numerical 
values  of  the  energy  and  momentum  relaxation  times  as  determined  in 
the  remainder  of  this  section.  These  processes  are  apparently  of 
minor  importance  in  comparison  to  those  in  the  first  two  collision 
categories  and  are  ignored  in  the  present  study,  although  their 
inclusion  offers  no  difficulty  in  principle. 
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2.2.3  The  Complete  Transport  Equations .  Adding  the  collision 


terms  whose  forms  have  now  been  determined  to  the  carrier,  velocity, 
and  energy  transport  equations  results  in 
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The  complete  energy  and  momentian  conserving  transport  model  consists  of 
seven  equations,  the  electron  and  hole  versions  of  Eqs.  2.33  through 
2.35  together  with  Poisson's  equation  for  the  electric  field, 

f  -  f  (»!-»•  P)  .  (2.36) 

where  is  the  net  density  of  ionized  impurities. 

2.2.4  Evaluations  of  Energy-Dependent  Parameters .  The 
collision  terms  in  Eqs.  2.33  through  2.35  contain  ionization  rates  and 
relaxation  times  which  are  as  yet  unknown  functions  of  carrier  energy. 
These  functions  must  be  given  numerical  values  before  the  energy  and 
momentum  conserving  transport  model  can  be  applied  to  device  simulation, 
although  the  form  of  the  model  is  independent  of  the  particular  pro¬ 
cedure  used  to  obtain  numerical  values.  The  method  of  evaluation 
>diich  will  be  used  here  will  be  to  determine  the  functions  using 
values  of  drift  velocity  and  ionization  rate  which  are  known 
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experimentally  as  functions  of  dc  electric  field.  A  mapping  of  the 
desired  quantities  onto  the  dc  field  is  obtained  by  writing  simplified 
versions  of  the  transport  equations  applicable  to  the  conditions  under 
which  experimental  measurements  are  made  and  substituting  in  the 
known  functions  of  field.  A  theoretically  determined  relationship 
between  energy  and  dc  field  is  used  to  complete  the  mapping  onto 
energy.  The  resulting  relationships  between  ionization  rates,  relaxa¬ 
tion  times,  and  energy  will  be  assumed  to  hold  under  all  conditions. 

This  procedure  for  evaluating  the  functions  of  energy  is  not  the 
only  possible  ^ne,  but  it  has  several  advantages.  It  makes  use  of 
comparatively  simple  equations  involving  experimentally  measured  values 
of  ionization  rate  and  drift  velocity,  and  guarantees  that  results  of 
the  model  will  agree  with  experiment  in  the  static  limit.  It  predicts 
forms  for  the  energy  and  momentum  relaxation  times  which  are  in  accor¬ 
dance  with  physical  expectations,  such  as  both  relaixation  times  be¬ 
coming  decreasing  functions  of  energy  in  the  energy  range  where 
impact  ionization  becomes  significant.  Finally,  although  the  procedure 
makes  use  of  certain  assumptions  about  the  forms  the  transport  equations 
can  take  in  the  presence  of  spatially  uniform  dc  electric  fields,  these 
are  all  confirmed  by  simulation  results  presented  in  Chapter  IV. 

2.2.I4.I  The  Static  Transport  Equations.  It  will  be  assumed 
that  the  circumstances  to  which  experimentally  measured  values  of 
electron  drift  velocity  and  ionization  rate  are  appropriate  include 
those  of  spatially  uniform,  dc  electric  field  and  negligible  hole 
concentration.  Under  these  conditions,  all  time  derivatives  and  hole 
ionization  terms  drop  out  of  the  transport  equations.  The  spatial 
derivatives  of  u  and  w  will  be  neglected.  (Simulation  results  for  the 


case  of  spatially  uniform,  dc  field  presented  in  Chapter  IV  support 
the  use  of  this  assumption.)  The  stated  conditions  of  zero  time 
variation  and  zero  spatial  variation  of  E,  u,  and  w  will  be  referred 
to  as  "uniform"  conditions,  under  which  the  transport  equations 
become 
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(2.39) 


The  two  terms  on  the  right-hand  side  of  Eq.  2.38  can  be  identi¬ 
fied  with  field-driven  drift  and  diffusion  down  the  carrier  concentration 
gradient.  Equating  the  first  of  these  with  the  conventional  drift 
velocity  yields 
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V 


(2.1*0) 


Since  v.  is  known  as  a  function  of  the  uniform  electric  field,  Eq.  2.1*0 
d 

relates  t  to  the  field  under  the  assumed  static  conditions. 

V 

The  right-hand  side  of  Eq.  2.37  expresses  the  same  quantity  as 
the  conventional  generation  rate  due  to  ionization  by  electrons. 
Equating  the  two  gives 

an  =  -  J  a«  ,  (2.1*1) 

q  n 

where  is  the  electron  current  density,  and  the  conventional  field- 
dependent  ionization  rate  has  been  labeled  with  an  asterisk.  Since 
u  is  the  average  electron  velocity  and  is  a  function  of  the  uniform 
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field,  Eq.  2. hi  is  equivalent  to 


a  =  uo* 


(2.i*2) 


This  maps  a  onto  E  under  uniform  conditons.  With  the  use  of  Eqs.  2.h0 
and  2.U2,  Eqs.  2.37  through  2.39  become 


(2.U3) 


^UT  - 

=  quET^  -  (w  -  -  mu^)a* 


(2.)+l*) 


Given  the  strength  of  the  electric  field,  v^  and  a*  are  known  in 

Eqs.  2.1*3  and  2.1*1*.  But  the  two  equations  are  still  in  terms  of  three 

unknowns,  u,  w,  and  t  ,  so  there  is  not  enough  information  to  map 

w 

T^,  T^,  and  a  onto  the  energy  w.  Equations  2.1*3  and  2.1*1*  with  any 
mapping  between  energy  and  field  will,  however,  guarantee  that  in  the 
limit  of  slow  changes  of  the  electric  field  in  time  and  space,  results 
from  the  energy  and  momentum  conserving  model  will  agree  with  those 
from  experiment.  The  next  step  is  to  find  a  physically  reasonable 
mapping  between  the  field  and  the  corresponding  mean  energy  for  momentum 
distributions  in  equlibriiim  with  the  field.  Several  calculations  of 
this  mapping  have  been  performed  in  the  course  of  development  of 
theories  for  the  field  dependence  of  the  ionization  rates . ® ^ ® ® 

2.2.1* .2  The  Uniform  Energy-Field  Relationship.  The  esti¬ 
mate  for  the  distribution  function  which  has  been  used  in  the  present 
work  was  developed  by  Wolff and  has 

N(c,E)  =  A  +  B  Ei(mc^/2F)  ,  0  <  c  <  c^ 

=  0  ,  c  >  ,  (2.U5) 
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vhere  A  and  B  are  normalizing  constants,  c  is  the  carrier  speed, 
is  the  speed  at  the  energy  threshold  for  impact  ionization,  and 
Ei  is  the  exponential  integral  function.  F  is  a  function  of  the  field 
defined  by 
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(qEX)^ 
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(2.1*6) 


where  X  is  the  mean  free  path  for  optical  phonon  emission  and  flin 
is  the  optical  phonon  energy.  X  and  the  ionization  threshold  energy 
are  chosen  so  that  the  ionization  rates  predicted  by  Wolff's  theory 
are  in  ’•easonable  agreement  with  measured  values.  Wolff  assumed  that 
ionization  takes  place  relatively  quickly  once  a  carrier  climbs  above 
the  energy  threshold,  so  that  negligibly  few  carriers  are  above  threshold 
at  any  one  time. 

Wolff's  distribution  cannot  apply  at  low  values  of  field  because 
it  makes  no  allowance  for  equilibration  between  the  carrier  and  lattice 
temperatures  as  the  field  approaches  zero.  It  will  be  assumed  that  the 
distribution  does  apply  when  the  field  is  at  least  300  kV/cm,  a  value 
at  which  the  static  ionization  rate  has  begun  to  be  appreciable.  In 
this  range,  the  desired  static  w-E  relation  is  just  the  second  moment 
of  Wolff's  distribution: 

c^ 

w(e)  =  f  ^  mc^N(c,E)c^  dc  ,  (2.U7) 

Jo 

where  the  differential  voltnne  is  a  spherical  shell  in  velocity  space. 

Ihe  integration  in  Eq.  2.47  must  be  performed  numerically  because  of 
the  presence  of  the  Ei  function  in  the  integrand. 

For  uniform  fields  of  less  than  300  kV/cm,  it  will  be  assiimed 
that  T  is  given  by  a  polynomial  in  the  field.  Knowledge  of  t  and 

V  w 
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3t  /3E  at  the  boundaries  of  the  interval  is  sufficient  to  determine 
w 

the  coefficients  of  a  third-order  polynomial,  t  and  3t  /3E  at 

w  w 

300  kV/cm  can  be  inferred  from  Eqs.  2.U3,  2.1+1*,  and  2.1+7.  The  re¬ 
maining  two  boundary  conditions  are  provided  by  examination  of  w  when 
the  field  is  close  to  zero.  Since  w  is  a  minimum  at  zero  field,  it 
must  be  a  function  of  even  powers  of  field  in  the  neighborhood  of 
zero.  It  will  be  assumed  that  the  carrier  distribution  is  not  signi¬ 
ficantly  heated  at  low  fields  so  that 


for  fields  close  to  zero.  Since  the  impact  ionization  rate  is  negli¬ 
gible  at  low  fields,  Eq.  2.1+1+  gives 


w  =  quEt  +  w  (2.1*9) 

w  o 

under  this  condition.  Combining  Eqs.  2.1+8  and  2.1+9  results  in 
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where  u  is  the  low-field  mobility.  The  fact  that  w  is  even  in  E, 
together  with  Eq.  2.1*8,  implies  that 


3t 

_ 5 

3E 


=  0 


E=0 


(2.51) 


Though  it  is  not  important  for  the  energy-field  relation,  it  is 
interesting  to  note  that  a  similar  analysis  involving  Eq.  2.1+3  leads 
to 
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This  can  also  be  seen  directly  from  Eq.  2.1+0. 
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With  T  (E)  determined  for  0  ^  E  <  300  kV/cm,  the  w-E  relationship 

V 

on  this  same  range  can  be  determined  from  Eqs.  2.1*3  and  2.1*1*.  It  would 
be  more  direct  to  assiime  that  w  itself  is  a  polynomial  in  E  on  this 
range,  but  this  could  result  in  w  not  monotonically  increasing  with  E, 
which  is  nonphysical. 

2. 2.1*. 3  The  Functions  of  Energy .  Given  the  full  range 
E-w  relationship,  the  relaxation  times  and  ionization  rate  can  be 
found  as  functions  of  w.  The  parameters  used  in  this  process  are 
shown  in  the  Appendix  for  both  electrons  and  holes.  They  include 
phonon  energy  and  creation  mean  free  path,  ionization  threshold  energy, 
low- field  mobility,  effective  mass,  and  parameters  for  the  phenomenolog¬ 
ical  relationships  between  drift  velocity,  conventional  ionization  rate, 
and  static  field.  The  lattice  temperature  is  taken  to  be  500°K.  The 
effective  mass  values  used  are  those  which  apply  at  the  energy  minima, 
though  they  could  be  changed  without  affecting  the  form  of  the 
transport  model.  Figure  2.1  shows  the  static  E-w  relationships  for 
electrons  and  holes,  while  Figs.  2.2  and  2.3  give  the  relaxation 
times  and  ionization  rates  as  functions  of  energy. 

2.3  Relationship  to  Other  Models 

In  Chapter  I,  it  was  mentioned  that  other  nonstatic  IMPATT 
models  have  been  attempted.  These  will  now  be  assessed.  The  condi¬ 
tions  under  which  the  energy  and  momentum  conserving  model  limits  to 
the  conventional  drift-diffusion  model  are  also  described. 

2.3.1  The  Energy  Conserving  Model.  Kafka  and  Hess'*^  developed 
an  IMPATT  model  which  includes  energy  conservation  and  energy-dependent 
ionization  rates,  together  with  conventional  field-dependent  veloci¬ 
ties.  Their  energy  transport  equation  is  similar  to  Eq.  2.35.  although 
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it  includes  a  nonzero  heat  flow  vector,  but  it  is  written  in  terms 
of  the  total  energy  of  electrons  and  holes  together,  so  that  the 
model  does  not  treat  the  electron  and  hole  energies  independently. 
Kafka  and  Hess  assert  that  their  model  becomes  equivalent  to  the 
conventional  one  when  the  flux  divergence,  pressure,  and  heat  flow 
terms  are  dropped  from  the  energy  transport  equation.  Since  they 
saw  significant  differences  in  their  simulation  results  depending 
upon  whether  these  terms  were  dropped,  they  concluded  that  the  terms 
represent  effects  not  allowed  for  by  the  conventional  model  which  are 
important  factors  in  the  operation  of  millimeter-wave  IMPATTs. 

Unfortunately,  the  reduced  form  of  the  Kafka  and  Hess  transport 
model  is  inconsistent.  The  pressure  term  in  the  energy  transport 
equation  does  not  drop  out  unless  the  carrier  concentration  gradient 
is  zero,  but  a  nonzero  gradient  must  in  fact  exist  when  impact  ioniza¬ 
tion  is  present.  Instead  of  disappearing,  the  pressure  term  in  the 
energy  transport  equation  of  the  reduced  model  should  take  the  form 
which  appears  in  Eq.  2.39* 

It  is  not  certain  that  this  oversight  affected  the  results 
obtained  by  Kafka  and  Hess,  but  it  may  have  done  so.  Under  static* 
conditions,  all  the  terms  in  the  energy  transport  equation  become 
small,  except  those  which  appear  in  Eq.  2.39.  These  remaining  terms 
determine  a  relationship  between  energy  and  field,  hence,  also 
betweeen  energy-dependent  ionization  rates  and  field.  In  the  reduced 


•The  term  "static"  Is  used  here  in  the  sense  as  defined  in  Chapter  I, 
where  it  refers  to  conditions  of  sufficiently  slow  space  and  time 
variation  of  the  electric  field  to  allow  the  carrier  momentum  dis¬ 
tribution  to  reach  equilibrium  with  the  field. 
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version  of  the  model,  one  of  the  terms  in  Eq.  2.39  is  missing,  and  a 
different  relationship  between  ionization  rate  and  field  results. 

What  Kafka  and  Hess  may  in  effect  have  done  was  to  perform  two 
essentially  static  simulations  with  different  sets  of  ionization  rates. 
This  possibility  is  consistent  with  the  facts  that  they  were  able  to 
make  their  results  from  the  two  models  agree  by  adjusting  the  ioniza¬ 
tion  rates,  and  that  their  plot  of  the  carrier  temperature  profile 
shows  relatively  constant  temperature  and  field  across  the  ionization 
region  of  their  device,  indicating  that  tlie  carrier  momentum  distri¬ 
bution  did  reach  equilibrium  with  the  field  in  precisely  the  region 
where  impact  ionization  was  taking  place. 

2.3.2  Energy  and  Momentum  Balance  Models.  Constant  and  co¬ 
workers  have  developed  two  Read-type  IMPATT  models  based  on  energy 

and  momentum  balance  relationships.  Their  first  model  applied  these 
relationships  to  the  IMPATT  drift  region  while  treating  the  ionization 
region  in  the  conventional  manner.^®  The  second  applied  the  energy 
balance  relationship,  with  energy-dependent  ionization  rates,  to  the 
ionization  region,  while  assuming  saturated  drift  throughout  the 
device. (This  was  done  partly  in  an  effort  to  confirm  the  effects 
of  avalanche  delay  due  to  energy  conservation  which  had  been  reported 
previously  in  connection  with  certain  results  from  the  present 
work.  )  Neither  of  the  energy  and  momentum  balance  models  treats 

the  entire  device  in  a  unified  way. 

The  balance  relationships  used  by  Constant  and  Salmer  are 
similar  to  those  proposed  by  Shur®’  and  are  equivalent  to 

^  (2.53) 
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where  is  assumed  to  be  constant.  These  relationships  are  Just 
approximate  versions  of  Eqs .  2.34  and  2.35-  They  omit  the  effects  of 
momentum  and  energy  flux  divergence,  of  the  pressure  terms,  and  of 
the  spatial  dependence  of  u,  w,  and  E.  The  balance  models  make  allow¬ 
ances  for  nonequilibrium  between  the  carrier  velocity  distribution  and 
the  electric  field  but  tend  to  overestimate  its  effects,  since  they 
assume  that  nonequilibrium  extends  iiniformly  over  the  entire  drift  or 
ionization  width,  ignoring  the  fact  that  carriers  will  tend  to  approach 
equilibrium  with  the  field  as  they  traverse  a  region  of  uniform  field. 
It  should  also  be  noted  that  diffusion  effects,  which  are  included  in 
even  the  conventional  static  model,  are  not  allowed  for  in  Eqs.  2.53 
and  2.54. 

2.3.3  The  Drift-Diffusion  Limit.  If  the  difference  between 
the  energy  and  momentum  conserving  model  and  the  conventional  drift- 
diffusion  model  lies  in  whether  nonequilibrium  conditions  are  allowed 
for,  then  the  two  models  should  be  equivalent  once  equilibrium  is 
reached.  This  is  indeed  the  case.  When  the  field  changes  slowly  in 
time  and  space,  the  velocity  transport  equation  reduces  to  Eq.  2.38. 
Substituting  Eq.  2.38  for  u  in  Eq.  2.33  gives 
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and  G  is  the  carrier  generation  rate.  Equation  2.55  is  just  the 

standard  drift-diffusion  equation  for  electrons  with  mobility  u  and 

diffusion  coefficient  D.  D  and  q  satisfy  the  Einstein  relation  in 

terms  of  T  .  D  can  be  calculated  as  a  function  of  field  from  Eq. 
c 

2.57  using  the  static  relationships  between  field,  energy,  and 
relaxation  time  found  previously.  Results  are  shown  in  Fig.  2.h.  This 
static  D(e)  relation  is  a  consequence  of  the  assumed  static  w(E);  one 
mapping  implies  the  other.  It  would  be  possible  to  assume  D(E)  and 
derive  w(e),  but  comparatively  little  i s  known  about  the  physical  situa¬ 
tion  in  Si  from  which  to  perform  ab  initio  construction  of  D(E), 
especially  at  high  field  strengths.  An  apparently  reasonable  choice  of 
D(E)  can  easily  lead  to  an  unreasonable  w{E),  such  as  one  where  w  is 
not  raonotonically  increasing  with  E  or  where  w  reaches  values  well  in 
excess  of  the  ionisation  threshold.  The  latter  would,  in  fact,  result 
from  the  D(E)  assumed  in  at  least  one  IMPATT  simulation  study  based  on 
the  drift-diffusion  model. 


2.U  Summary  and  Conclusions 

An  energy  and  momentum  conserving  transport  model  for  carriers 
in  Si  has  been  developed.  The  model  consists  of  transport  equations 
which  are  velocity  moments  of  the  phase-space  transport  equation.  The 
collision  terras  are  evaluated  by  requiring  that  under  conditions  of 
slow  tirae  and  space  variation  of  electric  field  the  results  of  the 
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FIG.  2.4  DIFFUSION  COEFFICIENT  VS.  ELECTRIC  FIELD  STRENGTH 

FOR  ELECTRONS  (SOLID  CURVE)  AND  HOLES  (DASHED  CURVE) 
IN  Si  AT  500°K. 


.  .r 


model  must  be  consistent  with  experimentally  measured  carrier 
velocities  and  ionization  rates. 

The  model  is  free  of  the  chief  limitation  of  the  conventional 
drift-diffusion  model  in  that  it  does  not  assume  equilibrium  between 
the  carrier  velocity  distribution  and  local  electric  field.  It  is 
also  more  general  and  self-consistent  than  all  previous  attempts  to 
produce  nonstatic  transport  models  applicable  to  IMPATTs. 


CHAPTER  III. 


NUMERICAL  METHODS  FOR  IMPATT  DIODE  SIMULATION 

This  chapter  describes  the  numerical  implementation  of  the 
ener/ry  and  momentum  conserving  transport  model  which  has  been  used  to 
simulate  the  operation  of  Si  IMPATT  diodes.  The  organization  of  the 
chapter  is  as  follows.  In  Section  3-1  a  set  of  normalizations  for 
finite-difference  approximations  to  the  transport  equations  is 
developed.  Section  3.2  considers  the  stability  and  accuracy  of  various 
finite-difference  forms  and  describes  the  form  which  has  been  chosen 
for  the  simulation  program.  Section  3.3  explains  how  the  program 
applies  spatial  and  temporal  boundary  conditions  to  the  finite- 
difference  equations,  and  Section  3.U  describes  the  findings  and 
conclusions  of  this  chapter. 

3.1  Normal i zat ions 

There  are  a  nximber  of  constants,  such  as  time  and  space  step 
length,  which  appear  repeatedly  in  finite-difference  approximations 
to  the  transport  equations  developed  in  Chapter  II.  if  the  finite- 
difference  equations  are  normalized  with  regard  to  these  constants, 
the  equations  are  simplified, and  efficiency  is  greatly  increased. 

This  section  describes  the  set  of  normalizations  which  has  been 
used  in  the  present  study. 

The  transport  equations  are  repeated  here  for  convenience: 

f  -  -If  .an.Rp  .  (3.1) 
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Their  peneral  finite-difference  approximations  are  given  by 
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where  6,  represents  any  finite-difference  in  time,  6  is  any  finite- 
difference  operator  in  space,  and  A  is  an  incremental  operator. 

(Only  constant  At  and  Ax  will  be  considered. )  When  the  same  notation 
is  used,  the  finite-difference  approximation  to  Poisson's  equation  is 


6  E 

X 

Ax 


^  (N,  -  N  +  p  -  n) 
e  d  a 


(3.T) 


In  Eq.  3.1*,  use  of  a  normalized  velocity 

At 


u  =  u 


AX 


(3.8) 
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and  normalized  generation  rates 


oAt ,6At 


(3.9) 


-k6~ 


so  that  it  is  convenient  to  define  a  normalization  for  carrier  and 
dopinp  concentration  by 

n  =  -a^n  .  (3.16) 

—  e 


Normalization  of  the  carrier  concentrations  has  no  effect  on  Eqs.  3.1^ 
through  3.6. 

Along  with  Eq.  3.10,  the  normalized  versions  of  Eqs.  3.^  through 

3.7  are 
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.  E  2 

-  u6  u  +  —  - 
—  X—  m  3n 


iii 

-  u6  w  +  — 
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’  (3.17) 

+  -  1 
n  ^ 

(3.18) 
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(3.19) 


Table  3.1  summarizes  the  normalizations.  Henceforth,  the  underscore 
notation  for  normalized  quantities  will  be  omitted,  and  when  it  is  not 
apparent  from  context  whether  the  quantities  referred  to  are  normal¬ 
ized  or  not,  this  will  be  stated  specifically. 


3.2  Finite-Difference  Operators 

There  are  many  possible  forms  for  Eqs.  3.10,  3.17,  and  3.18, 
corresponding  to  different  forms  of  the  operators  6^  and  6^,  and  to 
different  choices  of  time  levels  at  which  the  various  spatial  difference 
terms  are  evaluated.  An  optimum  finite-difference  form  will  be  one 
which  is  both  accurate  and  efficient,  but  the  goals  of  maximum 
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Table  3.1 


Normalizations  for  Finite-Difference  Forms 
of  the  Transport  Equations 


Cfuantity 


Average  velocity 
Average  energy 
Par-ticle  concentration 
Electric  field 
Relaxation  time 
Ionization  rate 


Symbol 

u 

w 

n 

E 

r 

ci,& 


Normalized  Quantity 
u  =:  u(At/Ax) 
w  =  (w/m){ At/Ax)^ 
n  =  n[ (q^ At^ ) /e] 
E  =  (qEAt^)/(Ax) 
T_  =  t/At 
Cl  *  aAt ,  RAt 


(  1 
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accuracy  and  maximtun  efficiency  are  incompatible.  This  section 
shows  how  a  compromise  between  the  two  has  been  reached. 

It  will  be  assumed  that  the  ratio  Ax/At  will  have  some  minimum 
value  which  will  be  given  by  the  form  of  the  Courant-Friedrichs-Lewy 
condition  which  applies  to  whatever  particular  numerical  form  is 
chosen  for  the  transport  equations.  If  this  ratio  is  always  given 
its  minimiim  value,  the  normalized  drift  and  pressure  terms  in  Eqs. 

3.10,  3.17  and  3.18  become  independent  of  Ax  and  At,  while  the  nor¬ 
malized  source  (carrier  generation  and  field)  and  relaxation  terms 
become  proportional  to  At. 

In  order  to  resolve  events  which  take  long  enough  in  relation 
to  the  RF  period  to  be  of  importance  to  device  operation.  At  itself 
must  be  sufficiently  small.  The  time  step  will  become  smaller  as 
frequency  increases,  and  vice  versa.  In  the  limit  of  very  small  At, 
the  drift  and  pressure  terms  will  dominate  the  normalized  finite- 
difference  equations,  and  in  the  limit  of  large  At,  the  field,  gener¬ 
ation,  and  relaxation  terms  will  dominate.  This  is  reflective  of  the 
shift  from  static  to  nonstatic  transport  as  device  speed  increases  and 
the  characteristic  time  scale  becomes  shorter.  In  a  comparatively  slow 
device,  static  equilibrium  between  the  field  and  relaxation  terms  ade¬ 
quately  describes  the  transport  of  carriers.  In  a  faster  device, 
the  nonstatic  effects  of  drift  and  pressure  on  the  momentum  and 
energy  become  important. 

It  is  convenient  to  begin  a  discussion  of  finite-difference 
approximations  to  Eqs.  3.1  through  3.3  at  the  limit  of  very  short 
time  scale  and  small  At,  when  the  normalized  source  and  relaxation 
terms  can  be  neglected.  If  it  is  to  be  useful  in  general,  any 
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approximation  must  be  useful  in  this  limit.  Therefore,  the  first 
part  of  this  section  is  concerned  with  the  sourceless  and  collision¬ 
less  version  of  the  finite-difference  model,  and  the  second  part  with 
the  further  considerations  which  arise  when  the  source  and  collision 
terras  are  restored  to  the  finite-difference  equations. 

3.2.1  Transport  in  the  Limit  of  Short  Time  Step.  The  source¬ 
less  and  collisionless  forms  of  Eqs.  3-10,  3.17,  and  3.18  are 


and 


6tW 


6^n  =  -  6^(nu)  , 

-  u6  u  -  -f-  6  nfw  -  ^  u^ 
X  3n  X  [_  '■  2 

-  u6  w  -  ■§-  S  nufw  -  4  u^ 
X  3n  X  r  2 


)] 

>]  • 


(3.20) 

(3.21) 


(3.22) 


Various  forms  of  the  difference  operators  will  now  be  considered. 

3. 2. 1.1  Forward- Time .  Upwind  Drift  Differences.  Using 
forward-time  differences  in  Eqs.  3.20  through  3.22  requires  less 
storage  than  centered  differences,  and  it  is  simplest  to  evaluate  all 
the  spatial  difference  terms  at  present  time.  Upwind  differencing 
tends  to  preserve  the  "transportive”*  property  of  drift ,  while 
centered  differencing  of  the  pressure  term  in  Eq.  3.21  will  reflect 
the  fact  that  acceleration  due  to  pressure  differences  acts  in  both 
the  upstream  and  downstream  directions.  The  pressure  terms  in  Eq. 
3.22  can  be  interpreted  as  keeping  track  of  the  work  done  by  carriers 
in  drifting  through  the  pressure  field,  so  upwind  differencing  is 
appropriate ®  A  reasonable  form  for  Eqs.  3.20  through  3.22  would 


•The  term  "transportive"  is  applied  here  as  it  is  defined  in 
Reference  1^9^  p.  6?. 
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therefore  appear  to  be 
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,  (3.25) 


where  subscripts  denote  position  in  the  space  mesh  and  superscripts 
position  in  time. 


3. 2. 1.2  Stability  Analysis.  The  stability  of  Eqs. 

3.23  through  3.25  can  be  examined  by  extending  the  usual  Fourier 
stability  analysis  to  three  variables  in  the  way  outlined  by 
Potter. The  solution  to  Eqs.  3.23  through  3.25  at  a  point  in  time 
is  a  vector  function  on  the  points  of  the  space  mesh.  This  function 
has  a  Fourier  decomposition  which  can  be  denoted  by 


M  _+  iJ3_ 
=  y  e  ® 

m=-M 


(3.26) 
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where 


6^  =  mir/M  ,  (3.2?) 

i  is  the  square  root  of  -1,  and  M  is  the  nvunber  of  points  in  the  mesh. 
The  numerical  error  ?  present  in  the  solution  has  a  similar  decom¬ 
position.  When  the  error  is  small,  it  changes  linearly  across  a 
time  step,  and  the  change  in  each  of  its  Fourier  components  can  be 
expressed  in  terms  of  an  amplification  matrix  [G]: 

^t+At  ^  ^  ])^t  _ 

m  mm 

The  identity  matrix  in  Eq.  3.28  accounts  for  the  effect  of  the 

forward  time  operator  used  in  each  of  Eqs.  3.23  through  3.25.  The 

particular  form  of  [G  ]  corresponding  to  these  equations  can  be 

™  _t 

found  by  substituting  the  Fourier  series  for  Cj  into  Eqs.  3.23  through 
3.25,  resulting  in  2M  +  1  component  equations  similar  to  Eq.  3.28. 
Linear  change  in  the  error  means  that  the  variable  coefficients  of 
the  spatial  differences  can  be  treated  as  constants  across  each  time 
step,  so  [G^]  is  given  by 

-iB  -iB 

-  u(l  -  e  ”*)  -  n(l  -  e  ”*)  0 

-  I  -  “(1  -  e  ")*|SSsInS,  -  |i  sin 

.f(v-|u")U  -e‘‘S  »(1  - 

(3.29) 
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stability  of  Eqs.  3.23  through  3.25  requires  that  no  component 
of  the  error  be  able  to  grow  in  magnitude  across  a  time  step.  Equa¬ 
tion  3.28  imnlies  that  this  will  be  true  if,  for  all  -  it  <  6  <  ii» 

m  ~ 

all  the  eigenvalues  X  of  [G  ]  satisfy 

m 


|1  +  X|*  ^  1  . 


(3.30) 


It  is  tedious  to  solve  the  characteristic  equation  of  [G  ] 

m 

directly,  but  there  is  a  simpler  way  of  getting  the  same  information. 

The  upper-left  entry  in  which  results  from  the  upwind  drift 

operator,  repeats  in  the  other  two  entries  on  the  diagonal,  so 

[G  ]  is  the  sum  of  two  simpler  matrices: 
m 


[G  ]  =  -  u(l-e  )[I]  +  [Gj  , 

m  m 


where 
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Any  eigenvector  of  [G^]  is  an  eigenvector  of  and  when  X'  is  an 

eigenvalue  of  [G'],  the  corresponding  eigenvalue  of  [G  ]  is  given  by 
m  IQ 


X  =  X’  -  u(l  -  e  )  . 


(3.33) 


-53- 


The  characteristic  equation  of  [G*]  is 

m 


x'[x'^  +  -|^X'{l-e  ®-i  sin  $  )  -  -^(w  -  ^  u^)(l  -  e  “)  sin  g  ] 

j  la  y  d  m 


(w  -  ^  u^)(l  -  e  ™)  sin  B  =  0  ,  (3.3*+) 

id  m 


which  results  in  X'  =0  and 
1 

X'  =  -  ^(1  -  3  -  i  sin  B)±4[(l-e  ™-i  sin  6 

2,33  m  3  m 

+  10i(w  -  i  u^)  sin  B„(l  -  e  .  (3.35) 

d  m 

X  corresponding  to  X'  satisfies  Eq.  3.30  for  any  value  of  B  > 
11  m 

X'  and  X'  are  somewhat  complicated  functions  of  B  «  hut  they  can  be 
2  3  m 

simplified  by  noting  that  w  is  nearly  always  large  in  comparison  to 

u^.  The  magnitudes  of  X’  and  X’  are  largest  when  6  =  v/2, 

2  3  HI 

when  the  two  eigenvalues  are  given  approximately  by 

.  -  |±  [i|v(i  -  1)]^  .  (3.36) 

When  the  plus  sign  is  used,  Eq.  3.30  becomes  approximately 

|l  -  1.33u  +  0.1+8/w  +  1.2i/w|2  <_  1  .  (3. 37) 

Equation  3.37  is  seldom  if  ever  satisfied  under  the  conditions 
which  occur  in  an  operating  IMPATT.  This  means  that  the  numerical 
method  given  by  Eqs.  3.23  through  3.25  is  unstable  in  the  limit  of 
small  At.  It  will  be  shown,  however,  that  when  At  is  sufficiently 
large,  the  method  can  be  stabilized  by  proper  treatment  of  the 
relaxation  terms,  so  that  it  has  been  possible  to  use  the  method  as 
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shown  in  Eqs.  3.23  through  3.25  in  much  of  the  diode  simulation  work 

performed  in  the  course  of  this  study. 

3.2.1. 3  Forward-Time ,  Centered  space  Differences .  In 

order  to  develop  a  stable  method,  it  is  useful  to  examine  another 

form  of  Eqs.  3.20  through  3.22.  The  reason  for  the  instability  of 

Eqs.  3.23  through  3.25  is  that  appears  in  the  real  part  of  X, 

where  it  adds  in  Eq.  3.37  to  the  number  one  which  appears  from  the 

forward-time  operator.  The  presence  of  the  quantity  1  -  e  in 

[C]  is  what  causes  this.  If  all  the  upwind  differences  which  contri- 
m  .  * 

-iBm 

bute  1  -  e  are  changed  to  centered  differences,  they  will  instead 

contribute  i  sin  3  ,  and  [O']  will  become 
m’  m 

= 

0  -  in  sin  B  0 

m 

2i/  1  „  2iu  .  „  2i  . 

-  -  ■^u'^)  sin  B  sin  B  -  sin  6 

3n  2  m  3  m  3  m 

2iu/  l2\.„  2i/  32\.„  2iu 

-  ^  )  sin  6„  -  )  sin  B  -  -7-  sin  B„ 

sn  d.  mq^  mom 

(3.38) 

The  eigenvalues  of  [G*]  are  now 

m 

X'  =0 
1 

and 

X^  ^  =  ±  i  sir.  .  (3.39) 

But  now  when  Bj^j  =  it/2,  Eq.  3.30  becomes 
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ll  ±  i[^(v  -  £  1  .  (3.40) 

Equation  3.40  can  never  be  satisfied,  so  the  modified  version  of  Eqs. 
3.23  through  3.25  is  no  more  stable  than  the  original. 

3. 2.1.4  The  Lax  Method.  There  is  a  stable  form  of 
Eqs.  3.20  through  3.22  which  is  explicit  and  involves  just  two  time 
levels.  It  is  based  on  a  nximerical  form  developed  by  Lax®®  for  the 
usual  equations  of  hydrodynamics,  which  are  similar  to  Eqs.  3.20 
through  3.22.  This  form  uses  centered  differencing  for  all  differ¬ 
ences  in  space  but  differs  from  the  one  just  described  in  that  it 
uses  a  modified  forward-time  operator: 

6  X  =  -  — (x^  +  x^  )  (3  4l) 

Vj  ^j  2^  j+i  j-i^  ’  u.'+i; 

where  x  is  any  of  the  normalized  transport  variables.  The  new  time 
operator  changes  Eq.  3.30  to 

Icos  6  +  XP  <  1  ,  (3.42) 

m  '  — 

and  centered  drift  changes  Eq.  3.33  to 

X  =  -  iu  sin  B  +  X'  .  (3.43) 

m 

For  the  Lax  method,  [G^]  and  its  eigenvalues  remain  as  in  Eqs.  3.38 
through  3.39,  so  Eq.  3.42  requires 

1^1  +  (^('^  -  <  1  .  (3.44) 

In  terms  of  unnormalized  quantities,  Eq.  3.44  requires 


(3.45) 
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Equation  3-^5  is  a  statement  of  the  Courant-Friedrichs-Lewy 
condition  as  it  applies  to  Eqs.  3.20  through  3.22.  The  right-hand 
side  is  the  sum  of  the  advection  speed  |u|  and  the  speed  at  which  a 
pressure  disturbance  can  propagate.®®  The  latter  is  analogous  to  the 
speed  of  sound  in  a  fluid.  This  speed  is  often  several  times  larger 
than  |u I  when  carriers  in  a  semiconductor  are  considered,  so  Eq.  3.^5 
restricts  At  to  be  several  tin^s  smaller  in  relation  to  Ax  than  does 
the  least  restrictive  condition  which  may  apply  to  simulations  using 
the  drift-diffusion  model,  i.e.,®® 


Ax  ^ 

At  -  ^d  ’ 


(3.U6) 


where  v,  is  the  static  drift  velocity.  (The  form  of  Eq.  3.^5  also 
d 

illustrates  the  usefulness  of  having  the  same  Ax /At  in  all  simula¬ 


tions,  since  the  stability  of  nxrmerical  methods  tends  to  depend  on 
this  ratio. ) 

The  stability  of  the  Lax  method  is  achieved  at  the  price  of  a 
comparatively  large  amount  of  numerical  diffusion.  The  method  intro¬ 
duces  a  spurious  diffusion  term  into  each  transport  equation,  with  a 
diffusion  coefficient  given  by"*® 


n  2At (  Ax^J 


(3.U7) 


where  u  is  the  unnormalized  drift  velocity.  Since  Eq.  3.^5  usually 
requires  that  Ax/At  be  at  least  5  x  10^  cm/s,  can  easily  exceed 
the  effective  diffusion  coefficients  derived  in  Chapter  II. 

3. 2. 1,3  Three- Level  Schemes .  Numerical  diffusion  can 
be  substantially  eliminated  by  the  use  of  a  three-time-level 
scheme.  A  number  of  such  schemes  have  been  developed  for  the 


> 


hydrodynamic  equations  and  are  generally  considered  to  be  variations 
on  a  scheme  originally  presented  by  Lax  and  Wendroff.®®  Each  of 
these  schemes  is  subject  to  the  same  stability  limitation,  as  given 
by  Eq.  3.^5,  as  the  Lax  method. The  original  Lax-Wendroff  scheme 
is  two  level,  but  is  somewhat  complicated.  Richtmyer®^  proposed 
an  equivalent  but  simpler  two-step,  three-level  scheme.  Richtmyer's 
method  gives  a  useable  solution  only  at  alternate  points  on  the  space 
mesh  and  only  on  alternate  time  steps,  so  a  variation  developed  by 
Burstein®®  has  been  chosen  for  use  in  the  simulation  program.  In 
this  method,  the  transport  equations  are  considered  as  a  single¬ 
vector  partial-differential  equation: 


i£  =  ii 

3t  ■  3x  ’ 


(3.U8) 


Equations  3.20  through  3.22  can  be  considered  to  be  in  the  form  of 
Eq.  3.^8  if  the  variable  coefficients  of  the  space  derivatives  in 
Eqs.  3.21  and  3.22  are  treated  as  constants  across  each  time  step. 

The  first  step  in  Burstein 's  method  uses  the  Lax  method  to 
find  a  trial  solution  at  the  half  time  and  half  space  step: 


;^+At/2  _  l/7t 

■  2^^J+i 


+  ft)  -  ^( 

j'  2Ax^ 


Tt 

J+l 


(3.it9) 


The  second  step  uses  these  centered  values  to  evaluate  t  and  advance 
the  solution  across  a  full  time  step: 


Burstein's  method  is  stable  provided  Eq.  3.^5  is  satisfied,  and  it 
introduces  little  spurious  diffusion.  Its  usefulness  is  obtained 
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at  the  cost  of  sacrificing  the  transportive  property  of  upwind 


drift  and  of  increased  computational  effort  in  comparison  to  a 
one-step  scheme.  Space  centered  differencing  does  give  some  ad¬ 
vantage  in  terms  of  program  simplicity  because  it  requires  no  special 
treatment  of  velocity  reversals. 

3. 2. 1.6  Other  Schemes.  It  should  be  noted  that  methods 
other  than  the  Lax  and  Lax-Wendroff  types  may  be  applicable  to 
Eqs.  3.20  through  3.22.  Richtmyer^®  describes  a  two  time  level  sbheme 
for  the  hydrodynamic  equations  which  advances  u  in  time  before  it 
does  n  and  w,  but  he  gives  an  apparently  erroneous  stability  analysis. 
It  is  unclear  whether  this  scheme  would  adapt  successfully  to  the 
transport  equations  used  in  this  work.  No  implicit  scheme  has 
achieved  any  significant  degree  of  acceptance  among  fluids  simu¬ 
lators,**®  though  a  workable  one  has  been  developed  by  Polezhaev.’® 

This  scheme  is  unfortunately  only  applicable  to  "supersonic"  flow, 
in  which  u  is  greater  than  the  square-root  quantity  in  Eq.  3.^5. 
Several  other  implicit  schemes  have  been  tried  on  the  computer  in 
the  course  of  the  present  work.  None  was  found  to  be  reasonably 
accurate  and  to  have  greater  stability  than  the  Burstein  method. 

3. 2. 2  Source  and  Relaxation  Terms. 

3. 2. 2.1  Carrier  Generation  and  Electric  Field  Terms.  The 
impact  ionization  term  represents  exponential  growth  of  the  carrier 
concentration,  with  growth  rate  a  or  B  (unnormalized).  In  the 
simulation  of  millimeter-wave  devices,  the  product  of  growth  rate 
and  time  step  is  much  less  than  one,  so  a  first-order  approximation 
of  the  exponential  growth  is  sufficiently  accurate: 

=  an^  +  Bp^  .  (3.51) 
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In  the  simulation  of  microwave  devices,  the  time  step  can  become 
large  enough  that  a  second-order  treatment  of  generation  is 
required,*®  but  this  will  ordinarily  occur  only  in  situations  where 
the  static  transport  model  is  adequate  for  describing  device 
behavior. 

The  terms  containing  the  electric  field  in  the  energy  and 
momentum  transport  equations  are  evaluated  at  present  time.  The 
coupling  of  the  field  terms  to  space-charge  density  through  Poisson's 
equation  introduces  restrictions  concerning  the  dielectric  relaxa¬ 
tion  time  and  the  Debye  length.  Preventing  numerical  overshoot  of 
the  charge  concentration  in  low-field  regions  requires  that  At  be 
shorter  than  the  dielectric  relaxation  time,*®  but  this  is  usually 
less  restrictive  than  Eq.  3.^5  when  Ax  is  chosen  for  reasonable 
spatial  resolution  in  a  millimeter-wave  device.  The  device  boun¬ 
daries  are  not  described  with  precision  when  Ax  is  greater  than  the 
Debye  length,  which  can  be  on  the  order  of  10“^  cm  in  contact 
regions,  but,  as  shown  by  results  given  in  Chapter  IV,  little  is 
gained  in  the  description  of  overall  device  behavior  by  making  Ax 
so  small. 

3. 2. 2. 2  Relaxation  Terms.  The  releucation  terms  represent 
exponential  decay  of  u  and  w.  Both  the  energy  and  momentum  equations 
have  the  form 

f  =  S-f  ,  (3.52) 

where  S  is  a  "source"  term  representing  the  influences  of  the  field 
aund  space  derivative  terms.  Tlie  exact  solution  of  Eq.  3-52  is 

f(t  +  At)  -  [f(t)  -  Sx]  +  ST  .  (3.53) 
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TVie  uimpler.t  finite-difference  approximation  to  Eq.  3.52  is  first 
order  in  time,  havinp  the  decay  term  e/aluated  at  time  t,  giving 

f^ 

^t+At  -  +  S^t  -  At  ,  (3.51*) 

Bosch  and  Thim*'  used  an  energy  and  momentum  conserving  model  to 
simulate  the  operation  of  transferred  electron  devices.  They  used 
a  second-order,  present-time  form  of  the  decay  term,  which  gives 

^t+At  ^  (ft  ^  ^  ^  ^  ^ 

2  2  T  2  T 


Another  possibility  is  to  perform  first-order  evaluation  of  the  decay 
term  at  advanced  time,  which  still  gives  an  explicit  solution.  The 
first-order,  advanced-time  approximation  to  Eq.  3.52  is 


ft+At 


ft+At 

f^  +  SAt - - -  At 

3  T 


f 3  +  SAt 
1  +  (At/r) 


(3.56) 


Setting  S  to  zero  in  Eqs.  3.53  through  3.56  gives  the  behavior 
in  each  of  any  ntimerical  error  which  may  be  present  at  the  beginning 
of  the  time  interval.  Error  is  plotted  in  Fig.  3.1  as  a  function  of 
At/x  for  each  of  the  forms  in  Eqs.  3.53  through  3.56.  The  figure 
shows  that  stability  of  f^  and  fj  requires  At  be  less  than  twice  the 
relaxation  time  because  the  magnitude  of  error  in  these  two  approxi¬ 
mations  grows  if  this  limit  is  exceeded.  This  would  be  a  severe 
restriction  on  At  in  actual  simulations,  since  the  momentum  relaxa¬ 
tion  time  in  Si  can  be  less  than  a  hundredth  of  a  picosecond.  No 
such  time-step  restriction  applies  to  f ^ ,  so  it  has  been  chosen  for 
use  in  the  diode  simulation  program,  f^  has  the  additional  advantage 
of  being  in  agreement  with  f(t)  as  given  by  Eq.  3.53  in  the  limit  of 

large  At.  This  is  not  true  of  f  or  f  . 

12 
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NORMALIZED  ERROR 


1.0 


At/r 

PIG.  3.1  NUMERICAL  ERROR  VS.  TIME  STEP  LENGTH  FOR  VARIOUS 

FORMS  OF  THE  DECAY  EQUATION.  [f(At);  EXACT; 

f  :  FIRST  ORDER,  PRESENT  TIME;  f  :  SECOND  ORDER, 
•  2 

PRESENT  TIME;  f  :  ADVANCED  TIME] 
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Advanced-time  relaxation  can  also  improve  the  stability 
which  is  associated  with  a  numerical  form  of  the  spatial  derivative 
terms.  This  can  be  seen  by  examining  its  effect  on  the  amplification 
matrix.  With  the  inclusion  of  source  terms  and  advanced-time  relaxa¬ 
tion  in  the  transport  equations,  Eq.  3.28  becomes 

=  [R]([I]  +  .  (3.57) 

where  R  is  a  relaxation  matrix  given  by 


[R] 


10  0 

0  [1  +  (1/t^)  +  6(p/n)]“^  0 

00  [1  +  (1/t  )  +  3(p/n)]"^ 

w 

(3.58) 


Since  the  diagonal  elements  of  [R]  are  less  than  or  equal  to  one,  the 
magnitudes  of  the  eigenvalues  of  the  right-hand  side  of  Eq.  3.58  are 
smaller  than  the  magnitudes  of  the  eigenvalues  of  [l]  +  [G^]  alone.  If 
the  latter  are  not  significantly  greater  than  one,  use  of  advanced¬ 
time  relaxation  can  stabilize  a  method  which  is  otheiT^ise  unstable. 

Equation  3.37  for  the  eigenvalues  associated  with  the  numerical 
method  given  by  Eqs.  3.23  through  3.25  indicates  that  keeping  the 
eigenvalue  magnitudes  close  to  one  requires  that  normalized  u  and  w 
be  sufficiently  small.  This  in  turn  requires  that  Ax/At  be  small,  so 
that  the  stability  requirement  for  the  method  in  Eqs.  3.23  through 
3.25  with  advanced-time  relaxation  amounts  to  a  Courant-Friedrichs- 
Lewy  condition.  In  practice,  it  has  been  found  that  the  method  is 
stable  if  Ax/At  is  greater  than  approximately  5  x  10^  cm/s  and  At  is 
greater  than  approximately  0.02  ps. 
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3. 2. 2, 3  Source  and  Relaxation  Terms  in  Burstein ' s  Method. 
If  requirements  for  spatial  resolution  dictate  that  Ax  be 
extremely  small,  the  CFL  condition  can  require  that  At  be  so  small  that 
the  normalized  diagonal  terras  in  the  matrix  [R]  in  Eq.  3.58  approach 
one.  If  this  happens,  advanced-time  relaxation  will  no  longer  provide 
numerical  stabilization,  and  a  Lax-Wendroff-type  scheme  for  the  spa¬ 
tial  derivative  terms  in  the  transport  equations  must  be  used.  Such 
schemes  advance  across  a  time  step  in  two  stages,  so  the  question 
arises  as  to  whether  the  source  and  relaxation  terms  should  appear  in 
both  stages  or  whether  the  effects  of  these  terms  over  a  time  step 
should  be  lumped  into  the  second  stage  only.  The  latter  approach  has 
been  followed  in  the  present  work.  Conceptually,  this  treats  the 
influences  of  the  pressure  and  drift  terms,  and  those  of  the  source  and 
relaxation  terms,  as  acting  in  parallel  across  a  time  step, just  as 
in  a  one-step  scheme.  The  approach  has  been  chosen  for  two 
reasons.  First,  the  time-centered  intermediate  result  obtained  in  the 
first  step  of  a  Leix-Wendroff  scheme  is  not  a  "true"  solution,  but 
only  an  estimate  upon  which  to  base  the  second  step.  Evaluation  of 
the  source  and  relaxation  terms  should  be  done  in  terms  of  the  gen¬ 
uine  solution,  which  is  not  available  as  a  time-centered  quantity. 
Second,  incorporation  of  source  and  relaxation  terms  only  in  the 
second  step  simplifies  the  method  somewhat. 

Burstein 's  method  with  source  and  relaxation  terms  in  the 
second  step  can  be  represented  as  follows.  The  first  step  remains 
exactly  as  in  Eq.  3.^9»  while  the  second  step  becomes 

f^^^^  =  +  +  ,  (3.59) 
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where  f  and  't'  are  defined  by  Eq.  3.*48,  [R]  is  the  relaxation  matrix 
Riven  by  Eq.  3.58,  and  S  is  a  vector  representing  the  contribution 
of  carrier  generation  and  the  electric  field. 

It  has  been  found  that  the  presence  of  source  terms  in  Burstein' 
method  can  cause  overshoot  in  the  solution  to  occur  near  inflow 
boundaries.  The  overshoot  can  be  eliminated  by  using  an  appropriate 
form  of  the  two-step  method.  The  sourceless  carrier  transport  equa¬ 
tion  can  be  written  in  two  numerical  forms  arising  out  of  the  compact 
and  expanded  forms  of  the  spatial  derivative  term: 


t+At  t  t,  t+At/2  t+At/2>  t,  t+At/2  t+At/2% 

n  -■  n  -  n  (u,,i  -  u,  i  )  -  u  (n  i  -  n  1  ) 

J  J  J  .J-5  j  J+5  j-5 


(3.62b) 


In  Eq.  3.21b  the  coefficients  n  and  u  of  the  space  derivatives  are 
centered  at  the  half  space  step. 

Equations  3.6la  and  3.6lb  are  numerically  equivalent,  but  there 
is  a  difference  between  Eqs.  3.62a  and  3.62b  which  gives  rise  to 
different  overshoot  properties.  At  the  inflow  boundary  in  an  IMPATT, 
carrier  concentration  and  velocity  can  both  change  rapidly  in  space 
because  carriers  enter  the  diode  from  a  low-field  region  where  they 
are  minority  carriers.  Once  inside,  they  accelerate  and  undergo 
impact  ionization,  and  the  resulting  changes  in  n  and  u  can  be  so 
rapid  as  to  cause  the  right-hand  side  of  Eq,  3.62a  to  be  negative  near 
the  boundary,  so  that  the  carrier  concentration  overshoots  past  zero. 
Equation  3.62b  contains  cross  products  between  n  and  u  at  different 
points  in  space,  so  the  tendency  to  overshoot  is  much  reduced. 


3.3  Initial  and  Boundary  Conditions 

The  transport  equations  determine  their  solution  to  within 
three  sets  of  conditions:  the  initial  conditions  on  the  simulation 
variables  throughout  the  space  mesh,  the  boundary  conditions  on  the 
simulation  variables  (usually  defined  at  the  edges  of  the  mesh ) ,  and 
the  relationship  between  terminal  current  and  voltage  which  is  deter¬ 
mined  by  the  interaction  between  the  diode  and  its  external  circuit. 
This  section  discusses  the  forms  and  methods  of  application  of  these 
conditions  in  the  simulation  program. 
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3.3.1  Initial  Conditions.  The  IMPATT  simulation  program  is 


"stand  alone"  in  the  sense  that  it  need  not  start  from  initial  condi¬ 
tions  provided  by  a  dc  solution  to  the  transport  equations.  Ibe 
program  can  either  set  up  initial  conditions  following  guidelines 
provided  by  the  user,  or  it  can  start  from  the  solution  obtained  at 
the  last  time  step  of  a  previous  simulation  run. 

The  use  of  an  arbitrary  starting  condition  which  might  never 
arise  in  the  course  of  actual  device  operation  might  appear  to  be 
in  contradiction  to  the  transport  equations,  but  any  starting  condi¬ 
tion  which  does  not  overspecify  the  initial  solution  is  mathematically 
possible.  Experience  with  the  simulation  program  has  never  turned  up 
a  situation  in  which  £in  arbitrary  set  of  starting  conditions  did  not 
evolve  rapidly  toward  a  physically  realistic  solution  as  the  simula¬ 
tion  progressed.  Care  must  only  be  taken  to  set  the  initial  conditions 
so  that  the  starting  transient  does  not  cause  unrealistically  large 
values  of  u  or  w  to  occur  momentarily,  since  this  can  violate  the 
condition  of  Eq.  3.^5. 

3.3.2  Spatial  Boundary  Conditions.  Even  though  the  transport 
equations  are  first  order,  centered  differencing  requires  that  boundary 
conditions  be  supplied  at  both  ends  of  the  space  mesh.  For  this 
purpose  a  psuedo  mesh  point  is  established  outside  each  end  of  the 
space  mesh,  and  it  is  at  these  two  points  that  the  boundary  conditions 
are  applied.  While  many  sets  of  boundary  conditions  are  possible, 
one  has  been  chosen  for  the  majority  of  the  present  work  which  com¬ 
bines  simplicity  with  a  reasonable  correspondence  to  the  limited 
knowledge  which  exists  concerning  the  conditions  in  the  contact 
regions  of  an  IMPATT  under  large-signal  operation.  At  an  inflow 
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boundary,  a  Dirichlet  condition  is  applied  to  each  unknown 
while  at  the  outflow,  a  Neumann  condition  with  constant  slope 
is  applied.  The  outflow  condition  is  chosen  to  maximize  the 
accuracy  of  the  finite-difference  scheme."*®  The  Dirichlet  inflow 
conditions  are 

n? 

\  •  <3-63) 

D 

=  0  (3.64) 

and 


(3.65) 


where  is  the  doping  concentration  in  the  contact  region.  The 
Inflow  conditions  describe  the  state  of  minority  carriers 
in  a  region  with  zero  electric  field.  They  might  seem  to  imply 
that  no  minority  current  can  enter  the  device,  but  it  is  possible 
for  a  finite  amount  of  minority  current  to  cross  the  device  boundary, 
which  is  located  halfway  between  the  pseudo  point  and  the  first 
actual  mesh  point.  This  ciirrent  is  small,  and  the  bulk  of  the 
reverse  saturation  current  is  provided  for  by  assuming  a  uniform 
rate  of  thermal  generation  throughout  the  device: 


G 


t 


(3.66) 


where  J  .  is  the  reverse  saturation  current  density,  and  L  is 
sat 

the  device  length.  Some  results  with  boundary  conditions  differ¬ 
ing  from  Eqs.  3.63  through  3.65  are  presented  in  Chai>ter  TV. 
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3.3.3  The  Device-Load  Interaction.  Updating  the  electric 
field  at  each  time  step  requires  that  some  relation  between  the 
terminal  voltage  and  current  be  supplied.  One  way  of  doing  this  is 
to  convolve  the  past  values  of  current  with  the  impulse  response  of 
an  external  circuit  after  every  time  step.^^  This  has  the  obvious 
drawback  of  requiring  that  many  convolution  integrals  be  evaluated 
in  the  course  of  a  simulation.  A  simple  way  of  centering  the  current 
in  the  upcoming  time  step  using  present  terminal  voltage  exists,^® 
but  it  does  not  allow  the  use  of  simple  capacitive  branches  (such  as 
an  RF  source  with  dc  block)  in  the  external  circuit.  Bauhahn 
and  Haddad®"*  imposed  a  sinusoidal  RF  voltage  at  the  device 
terminals.  While  their  method  requires  no  circuit  tuning  or  source 
adjustment  in  order  to  obtain  a  desired  RF  terminal  voltage  amplitude, 
it  does  require  iteration  on  the  dc  voltage  in  order  to  obtain  a 
desired  dc  current,  and  it  does  not  allow  the  simulation  of  genuine 
transients  in  connection  with  realistic  external  circuits.  A  state 
space  approach^ ^  to  the  device-circuit  interaction  has  been  found  to 
be  the  most  useful  in  the  present  work. 

The  load  routine  used  in  the  simulation  program  is  designed 
to  impose  a  sinusoidal  RF  voltage  on  the  device  terminals  in  a 
self-consistent  manner.  The  circuit  model  used  by  the  routine  is 
shown  in  Fig.  3.2.  The  diode  is  represented  by  a  particle  current 
source  in  parallel  with  the  diode  depletion  capacitance 
(normalized  to  unit  area).  The  external  circuit  consists  of  a  dc 
current  source  in  parallel  with  an  RF  voltage  source  series 

resistance  R,  and  blocking  capacitor  C.  The  strength  of  at  any 
point  in  time  can  be  found  from 
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q(u^n  +  Upp)  dx 


(3. 67) 


where  L  is  the  device  length. 

The  time  evolution  of  the  blocking  capacitor  voltage  and 
the  diode  voltage  are  given  by 


and 
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(3.69) 


(3. TO) 


Integrating  Eq.  3.68  and  3.69  across  a  time  step  using  the  trapezoidal 
rule  gives  an  equation  for  and  at  the  end  of  the  time  step: 
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In  Eq.  3.71»  is  known  from  Eq.  3.67. 
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If  the  external  series  resistance  R  is  chosen  to  be  zero,  the 


time  evolution  of  the  diode  voltage  is  given  by  Eq.  3-69,  with 


J 


(3.72) 


With  the  use  of  trapezoidal  rule,  Eqs.  3-69  and  3-72  give 


j'  -  .  C(v‘;«  .  v'p)l/(C  .  c^) 


(3.73) 


when  the  series  resistance  is  zero. 

The  form  of  the  circuit  in  Fig.  3.2  is  such  as  to  suppress 
unwanted  harmonic  components  of  the  terminal  voltage.  Bias  oscilla¬ 
tions’®  and  subharmonic  instabilities”  can  usually  be  prevented  by 
choosing  an  appropriate  ratio  between  C  and  C^.  The  series  resistance 
serves  to  damp  out  transients  more  quickly  them  they  would  otherwise 
decay.  Usually,  a  particular  steady-state  RP  voltage  amplitude  is 
desired  at  the  device  terminals.  The  terminal  amplitude  will  differ 
somewhat  from  that  of  the  RF  voltage  source,  but  can  be  bro\ight  to 
the  desired  value  by  adjustment  of  the  source  amplitude. 


3.U  Conclusions 

Development  of  an  IMPATT  simulation  program  based  on  finite- 
difference  approximations  to  the  energy  and  momentum  conserving 
transport  equations  requires  careful  choice  of  the  numerical  methods 
used.  The  methods  described  in  this  chapter  are  efficient  and  accur¬ 
ate  solutions  to  the  numerical  problems  associated  with  various  parts 
of  the  transport  equations. 

Explicit  finite-difference  schemes  incorporating  forward¬ 
time  differences  are  likely  to  be  unstable  in  the  limit  of  very 
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sliort  time  steps.  The  Lax  method  has  a  well-defined  stability 
requirement  which  is  acceptable  from  the  point  of  view  of  efficiency, 
but  the  method  introduces  an  undesirably  large  amount  of  numerical 
diffusion.  Ibis  diffusion  can  be  substantially  eliminated,  without 
sacrificing  stability,  by  the  use  of  a  scheme  of  the  Lax-Wendroff 
type. 

Both  the  first-  and  second-order  present-time  forms  of  the 
relaxation  terms  impose  a  severe  stability  restriction  on  the  time 
step.  The  first-order  advanced-time  form  imposes  no  such  restriction, 
and  its  asymptotic  behavior  is  the  same  as  that  of  the  exponential 
decay  process  which  has  been  used  to  model  the  effects  of  collisions. 

The  simulation  program  makes  use  of  initial  and  boundary 
conditions  which  are  simple,  stable,  and  consistent  with  reasonable 
assumptions  about  the  conditions  which  occur  in  an  actual  device. 
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CHAPTER  IV.  SIMULATION  OF 


MILLIMETER-WAVE  IMPATT  OPERATION 

This  chapter  presents  results  of  several  series  of 
investigations  which  have  been  performed  using  the  simulation 
program  described  in  the  preceeding  chapter.  The  first  section 
presents  dc  and  large-signal  results  for  situations  of  electric 
field  with  very  small  spatial  variation.  These  results  permit  clear 
identification  of  certain  overshoot  and  relaxation  phenomena  which 
are  present,  but  less  clearly  identifiable,  in  simulations  of  more 
realistic  IMPATT  structures.  The  results  also  demonstrate  the 
appropriateness  of  certain  assumptions  made  in  Chapter  II.  In 
Section  k,2  results  obtained  using  the  energy  and  momentum  con¬ 
serving  simulation  are  compared  with  results  obtained  using 
conventional  drift-diffusion  simulation  for  the  same  situations. 
Systematic  differences  are  observed  and  discussed. 

Sections  U.3  and  U.1+  are  concerned  with  numerical  experi¬ 
ments  to  establish  the  importance  of  particular  physical  mechanisms 
to  overall  device  behavior.  Section  4.3  is  concerned  with  carrier- 
inflow  boundary  conditions  and  considers  the  effects  of  highly 
doped  contact  regions  and  injection  of  "hot"  carriers  at  boundaries. 
Section  4.4  examines  the  effects  of  the  cooling  of  each  carrier 
distribution  due  to  impact  ionization  Initiated  by  the  other 
carrier  type. 

Section  4.5  is  concerned  with  the  implications  of  the  present 
study  for  the  potential  performance  of  millimeter-wave  Si  IMPATTs. 


Performance  limitations  due  to  nonstatic  carrier  transport  phenomena 
and  to  parasitics  external  to  the  active  diode  are  discussed.  A 
criterion  for  estimating  the  upper  frequency  limit  for  useful 
operation  of  the  IMPATT  mode  in  any  material  is  explained;  this 
upper  limit  is  estimated  to  be  approximately  500  GHz  for  Si  IMPATTs. 

4.1  "Flat  Field”  Results 

Simulation  of  pn  Junction  devices  with  spatially  uniform 
electric  field  in  avalanche  breakdown  is  useful  for  gaining  an  under¬ 
standing  of  the  behavior  of  the  solution  of  the  energy  and  momentum 
transport  equations  in  the  presence  of  time  and  space  variations  of 
the  electric  field.  With  flat  fields  a  dc  solution  shows  the  behavior 
of  the  transport  quantities  near  a  spatial  field  step  (present  at 
each  spatial  boundary)  in  the  absence  of  temporal  variations.  A 
large-signal  solution  shows  the  behavior,  away  from  the  sp  •‘■ial 
boundaries,  of  the  transport  quantities  under  time-varying,  spatially 
uniform  conditions.  The  responses  to  spatial  and  temporal  variation 
of  the  field  can  thus  be  observed  separately. 

Simulation  results  have  been  obtained  for  flat- field  situations 
in  devices  land  0.3  pm  long,  the  doping  profiles  of  which  are  shown 
in  Fig.  4.1.  The  doping  densities  are  chosen  to  approximately  com¬ 
pensate  the  space  charge  of  mobile  carriers  at  a  dc  bias  current 
density  of  10**  A/cm^.  Figures  4.2  emd  4.3  show  dc  carrier  concen¬ 
tration,  current  density,  average  velocity,  average  energy,  and 
electric  field  profiles  for  the  two  devices.  The  conventions  used 
in  the  figures  are  used  in  all  results  presented  in  this  chapter: 
positive  electric  field  and  hole  velocity  are  taken  from  right  to 
left,  and  positive  electron  velocity  from  left  to  right.  Electrons 
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FIG.  4.2  ELECTRIC  FIELD  (SOLID  LINES)  AND  CONCENTRATION, 


CURRENT  DENSITY,  AVERAGE  ENERGY,  AND  AVERAGE 
VELOCITY  FOR  ELECTRONS  (MINUS  SIGNS)  AND  HOLES 
(PLUS  SIGNS)  IN  THE  1-Min,  FLAT  FIELD  DIODE  UNDER 
DC  CONDITIONS.  (J^  =  10"  A/cm*) 
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FIG.  It. 3  ELECTRIC  FIELD  (SOLID  LINES)  AND  CONCENTRATION, 
CURRENT  DENSITY,  AVERAGE  ENERGY,  AND  AVERAGE 
VELOCITY  FOR  ELECTRONS  (MINUS  SIGNS)  AND  HOLES 
(PLUS  SIGNS)  IN  THE  0.3-pm,  FLAT  FIELD  DIODE 


UNDER  DC  CONDITIONS.  (J^^  =  10"  A/cm*) 
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enter  each  device  from  the  left-hand  boundary  with  concentration, 
average  velocity,  and  average  energy  as  described  in  Chapter  III; 
holes  enter  from  the  right. 

The  e' ectron  and  hole  concentrations  in  the  figures  increase 
monotonically  as  the  carriers  traverse  each  device.  This  is  the 
usual  situation  in  a  pn  Junction  under  reverse  breakdown.  The 
velocity  profiles  show  that  the  velocities  overshoot  their  conven¬ 
tional  "saturated"  values  as  carriers  enter  each  device,  approaching 
the  conventional  values  only  after  the  carriers  have  traveled  some 
distance  past  the  inflow  boundary.  These  velocity  overshoots  are 
opposed  rather  than  assisted  by  diffusion  down  the  carrier  concen¬ 
tration  gradients.  Carrier  energy  is  low  at  inflow  boundaries  and 
substantially  "catches  up"  to  the  field  over  the  velocity  overshoot 
distance.  Once  they  have  reached  equilibrium  with  the  field,  velocity 
and  energy  do  not  change  through  the  remainder  of  the  device  length, 
even  though  there  is  a  concentration  gradient  present.  This  is  in 
accordance  with  the  assinnption  of  "uniform"  conditions  which  was 
made  in  Chapter  II. 

Velocity  overshoots  similar  to  those  shown  in  Figs.  U.2  and 
lt.3  will  be  present  in  all  of  the  energy  and  momentum  conserving 
results  presented  in  this  chapter.  The  overshoots  occur  when 
carriers  enter  the  depletion  region  of  a  reverse-biased  diode.  In 
doing  so,  they  pass  throiigh  an  abrupt  step  in  field  strength,  since 
the  Inflow  boundary  conditions  approximate  conditions  in  a  low-field 
region.  Immediately  downstream  of  the  step  in  field,  the  carriers 
have  lower  average  energy  than  they  would  have  in  equilibrium  with  the 
field.  At  lower  energy  the  electron  and  hole  relaxation  times  are 


A0-A13a  901 

UNCLASSIFIED 


COMPUTf*  MOOIklNO  Of  «ILI IMiT|»-VAVI  IMPATT  DIODES  OME 
OF  A  SINUS  OF  NfP.  lOl  MICNIOAN  UNIV  ANN  AllSoN 
CLECTNON  PHYSICS  LA«  A  «  FMtLICH  NOV  ta  tirTsa 
AFNAl-TR-aa- 1107  f  SM  tft-DI  14)9  f/O  »/a 


microcopy  resolution  test  chart 


longer,  so  the  releucation  terms  in  the  velocity  transport  equations 
become  smaller,  and  velocity  overshoot  results. 

The  figures  suggest  that  the  carrier  velocities  respond  much 
more  quickly  to  changes  in  the  field  them  do  the  carrier  energies. 

Each  energy  can  be  seen  to  require  some  distance  over  which  to 
build  up  in  response  to  the  inflow  field  step.  In  contrast,  the 
corresponding  velocity  overshoots  almost  immediately,  and  the  distance 
over  which  the  overshoot  extends  corresponds  to  that  over  which  the 
energy  is  increasing.  This  indicates  that  velocity  transients  are 
due  primarily  to  the  energy  dependence  of  the  velocity  relaxation 
times.  The  situation  of  low  energy  and  velocity  overshoot  which 
occurs  at  each  inflow  boundary  in  Figs.  U. 2  and  U.3  will  be  referred 
to  as  spatial  lag,  i.e.,  the  lagging  of  energy  behind  a  spatial  var¬ 
iation  in  the  electric  field. 

Figures  U.U  and  U.5  show  profiles  of  electric  field,  average 
energy,  and  average  velocity  in  the  l-ym  flat-field  device  at  various 
points  in  a  large-signal  BF  cycle.  The  RF  terminal  voltage  is  sinu¬ 
soidal  with  a  frequency  of  300  GHz  and  an  amplitude  of  15  V,  which  is 
Just  under  half  the  dc  bias  voltage  of  32  V.  The  figures  show  condi¬ 
tions  inside  the  device  at  0-,  90- ,  l80-  and  270-degree  phase  in 
the  RF  cycle,  the  phase  angles  at  which  the  RF  terminal  voltage  or 
its  time  derivative  reaches  an  extremum.  These  figures  show  evidence 
of  spatial  lag  near  the  contacts.  Just  as  did  Figs,  k.2  and  U.3,  but 
now  the  region  of  interest  is  the  center  of  the  device,  where  the 
energies  and  velocities  can  be  seen  to  be  spatially  constant. 

The  carrier  energies  in  Fig.  U.U  rise  and  fall  under  the 
influence  of  the  electric  field.  Figure  U.5  shows  that  the  velocities 
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FIG.  4.4  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  ENERGY 
(MINUS  SIGNS)  AND  HOLE  ENERGY  (PLUS  SIGNS)  IN 


THE  1-vm,  FLAT  FIELD  DEVICE  AT  VARIOUS  POINTS 
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FIG.  4.5  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  VELOCITY 


(MINUS  SIGNS)  AND  HOLE  VELOCITY  (PLUS  SIGNS)  IN 


THE  1-nm,  PLAT  FIELD  DEVICE  AT  VARIOUS  POINTS  IN 
AN  RF  CYCLE.  (V^^  =  15  V  AND  =  10"  A/cm*) 


-82- 


depart  from  their  saturated  values  of  10^  cm/s  and  8,5  x  10®  cm/s 
for  holes  and  electrons,  respectively.  The  departures  occur  because 
of  time  lag  between  the  local  carrier  energy  and  electric  field. 

At  the  beginning  of  the  cycle,  field  strength  is  increasing  with 
time,  so  as  energies  lag  behind  the  field  the  velocities  overshoot 
for  the  reasons  discussed  in  connection  with  Figs,  k.2  and  At 

l80  degrees,  an  inverse  process  takes  place.  The  field  strength  at 
this  point  is  decreasing  with  time,  so  the  carrier  energies,  which 
lag  behind,  are  larger  than  those  which  would  occur  in  carrier-field 
equilibrium.  At  larger  energies  the  momentum  relajcation  times  are 
shorter,  so  the  velocities  undershoot  their  conventional  saturated 
values.  At  90  and  270  degrees  the  velocities  return  to  their 
satiurated  values  because  the  energies  substantially  "catch  up"  to  the 
slowly  changing  field  strength. 

Further  evidence  for  time  lag  between  field  strength  and 
carrier  energy  is  given  by  Fig.  U.6.  This  shows  the  carrier  ener¬ 
gies,  carrier  velocities,  and  field  strength  at  the  center  of  the 
1-um  flat- field  device  as  fxuictions  of  phase  over  one  RF  cycle. 

Time  lag  between  the  extrema  of  field  and  energy  is  apparent  in  the 
figure,  as  are  the  velocity  overshoots  and  undershoots  which  occur 
in  times  of  increasing  or  decreasing  field  strength.  The  time 
between  an  extremum  in  field  and  the  corresponding  one  in  energy 
is  approximately  1  ps.  The  lag  between  carrier  energy  and 
temporal  variation  of  the  electric  field  will  be  referred  to  as 


temporal  lag. 
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FIG.  4.6  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  ENERGY  (MINUS  SIGNS),  ELECTRON  VELOCITY  (MINUS 
SIGNS),  HOLE  ENERGY  (PLUS  SIGNS)  AND  HOLE  VELOCITY  (PLUS  SIGNS)  AT  THE  CENTER  OF  THE 
1-um,  FLAT  FIELD  DEVICE  VS.  RF  PHASE  ANGLE.  (V  =  15  V  AND  J  =  lO-  A/cm") 


h.2  Comparisons  Between  Conventional  and  Energy  and  Momentum 
Conserving  Results 

The  results  described  in  the  preceeding  section  indicate 
that  spatial  and  temporal  lag  between  carrier  energy  and  electric 
field  strength  occur  in  results  from  simulation  based  on  the 
energy  and  momentiim  conserving  transport  model.  This  lag  is  not 
accounted  for  in  conventional  simulation  based  on  drift  and  diffu¬ 
sion,  so  the  nature  and  degree  of  its  effect  on  IMPATT  operation  can 
be  determined  by  comparing  results  from  conventional  simulation*"  with 
those  from  energy  and  momentum  conserving  simulation.  Comparisons 
have  been  made  for  three  double-drift  IMPATT  structures  with  lengths 
of  1,  0.5,  and  0.3  pm.  The  three  doping  profiles  are  shown  in 
Fig.  it. 7.  They  were  chosen  so  that  the  devices  will  be  strongly 
punched  through  at  breakdown,  and  so  that  the  electron  and  hole 
drift  transit  times  will  be  approximately  equal  in  each  device.  The 
conventional  drift  velocities,  ionization  rates,  and  diffusion 
coefficients  used  are  given  in  Appendix  A  e^’d  Fig.  2.k. 

Figures  b,8,  it. 9,  and  U.IO  show  large-signal  admittance  restilts 
from  the  two  types  of  simulation  for  the  1-,  0.5-  and  0.3-pm  struc¬ 
tures,  respectively.  Table  it.l  lists  a  number  of  representative 
operating  characteristics.  There  is  very  little  difference  between 
the  sets  of  data  for  the  1-ym  structure.  A  consistent  difference 
between  corresponding  admittance  data  points  can  be  seen  in  the 
case  of  the  0.5-pm  structure}  and  in  the  case  of  the  shortest 
structure,  this  difference  is  more  pronounced.  Apparently  the  extra 
physical  effects  allowed  for  in  the  energy  and  momentum  conserving 
transport  model  have  a  significant  effect  on  the  operation  of 
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FIG.  U.7  DOPING  PROFILES  FOR  (a)  l-ym,  (b)  0.5-um 
AND  (c)  0.3-mn  DOUBLE-DRIFT  IMPATTS. 
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FIG.  U.8  ENERGY  AND  MOMENTUM  CONSERVING  (SOLID  LINE) 

AND  CONVENTIONAL  (DASHED  LINE)  G-B  RESULTS  FOR 

THE  1-um  DOUBLE-DRIFT  DEVICE.  (V„„  =  10  V  AND 

RF 

^  A/cm  *) 
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FIG.  4.9  ENERGY  AND  MOMEJITUM  CONSERVING  (SOLID  LINE) 
AND  CONVENTIONAL  (DASHED  LINE)  G-B  RESULTS 
FOR  THE  0.5-W  DOUBLE-DRIFT  DEVICE. 

(Vj^  =  8  V  AND  =  10*  A/ cm*) 
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FIG.  U.IO  ENERGY  AND  MOMENTUM  CONSERVING  (SOLID  LINE) 


AND  CONVENTIONAL  (DASHED  LINE)  G-B  RESULTS  FOR 
THE  0.3-Mm  DOUBLE-DRIFT  DEVICE.  (Vj^  =  6  V  AND 
J.  =  1.5  X  10*  A/cm*) 


Table  U.l 


Sample  Results  from  Large-Signal  Simulation 


Device 

(urn) 

•'dc 

(A/cm^  ) 

Frequency 

(GHz) 

V 

RF 

ill 

1 

6  X  10** 

•  100 

10 

0.5 

1  X  10® 

l4o 

8 

0.3 

1.5  X  10® 

200 

6 

G 

(mhos /cm* ) 

n 

(percent) 

V 

-  2.9  X  10® 

8.2 

30 

(-  2.7  X  10®)» 

(7.8)« 

(29.9) 

-  6.5  X  10® 

10. U 

20.3 

S  5.3  X  10®) 

(Q.k) 

(19.9) 

1.1  X  lO** 

8.6 

15.6 

(8.1  X  10®) 

(6.3) 

(15.0) 

^Results  from  conventional  simulations  are  shown  in  parenthesis. 


double-drift  Si  IMPATTs  for  device  lengths  of  approximately 
0.5  um  or  less. 

The  results  in  Table  U,1  show  that  the  energy  and  momentum 
conserving  simulation  generally  predicts  larger  dc  voltage  at  a 
given  operating  point  than  does  the  conventional  simulation.  This 
is  probably  an  effect  of  spatial  lag.  In  the  energy  and  momentum 
conserving  simulation,  the  low  energies  at  inflow  boundaries  result 
in  lower  rates  of  impact  ionization  than  would  be  predicted  by  the 
conventional  transport  model.  Lower  rates  near  the  boundaries  must 
be  compensated  for  by  larger  rates  in  the  center  of  a  device,  re¬ 
quiring  an  increase  in  field  strength  and  dc  terminal  voltage.  The 
occurrence  of  spatial  lag  in  a  double-drift  device  is  shown  by 
Fig.  4.11,  which  gives  energy, velocity,  and  field  profiles  in 
the  0.3-ym  double-drift  device  under  dc  conditions.  Carriers  enter¬ 
ing  the  device  can  be  seen  to  undergo  velocity  overshoot. 

The  way  in  which  the  admittance  results  in  Figs.  4.9  and 
4.10  differ  is  initially  surprising  if  one  intuitively  expects  that 
the  Inclusion  of  higher  order  transport  effects  in  an  IMPATT  model 
will  result  in  poorer  device  performance.  Where  the  results  differ 
significantly,  the  energy  and  momentum  conserving  simulation  con¬ 
sistently  predicts  larger  negative  conductance,  though  somewhat  lower 
optimum  frequency,  than  does  the  conventional  simulation.  There  are 
several  mechanisms  which  may  tend  to  increase  negative  conductance. 
One  of  these  is  delay  in  the  impact  ionization  process  due  to  energy 
lag.  This  would  cause  carrier  injection  to  occur  later  in  the  RF 
cycle.  Another  mechanism  is  velocity  overshoot  and  undershoot  on 
the  part  of  drifting  carriers.  Iliis  might  contribute  to  negative 
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conductance  by  giving  the  induced  current  waveform  a  more  favorable 
shape,  provided  overshoot  and  undershoot  were  to  occur  at  the  proper 
points  in  the  cycle.  Changes  in  velocity  might  also  affect  the  eunomt 
of  time  that  carriers  spend  in  the  ionization  region, thereby  affect¬ 
ing  the  shape  of  the  injected  pulse. 

Figures  U.12  through  U.l6  provide  insight  into  the  mechanisms 
which  cause  increased  negative  conductance  in  the  energy  and  momentum 
conserving  results.  The  figures  present  simulation  results  from  the 
0.3-um  device  operating  at  a  frequency  of  300  GHz  with  an  RF  amplitude 
of  10  V.  Figure  U.12  is  a  set  of  plots  of  injected  and  induced 
current  waveforms  resulting  frcan  the  two  types  of  simulation.  The 
sinusoidal  terminal  voltage  is  shown  for  phase  reference.  The 
injected  current  depicted  in  the  figure  is  the  integral  over  the 
device  length  of  the  instantaneous  electron  and  hole  generation  rates, 
taken  at  each  time  step.  The  induced  current  is  the  component  of 
terminal  current  which  flows  because  of  carrier  motion  in  the  interior 
of  the  device. 

Figure  U.12  shows  that  the  injected  current  waveforms  are  very 
similar  in  shape,  but  the  injected  current  in  the  energy  and 
momentxim  conserving  result  is  delayed  by  approximately  10  degrees 
(or  0.1  ps)  in  comparison  to  the  conventional  result.  A  similar 
delay  also  appears  in  the  plot  of  the  terminal  currents ,  along  with 
some  difference  in  shape.  Injection  delay  tends  to  increase 
negative  conductance  and,  by  shortening  the  optimum  transit  angle, 
tends  to  lower  the  optimum  IMPATT  operating  frequency.  Fig¬ 
ures  U.13  and  U.l6  show  profiles  of  carrier  concentration,  electric 
field,  average  velocity  and  average  energy  at  fo\ir  points  in  the  RF 
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FIG.  U.12  TERMINAL  VOLTAGE,  INJECTED  CURRENT,  AND 
INDUCED  CURRENT  WAVEFORMS  FOR  0.3-wm 
DOUBLE-DRIFT  IMPATT.  (a)  ENERGY  AND 
MOMENTUM  CONSERVING  SIMULATION  AND 
(b)  CONVENTIONAL  SIMULATION,  (f  =  300  GHz, 
Vpp  =  10  V  AND  =  1.5  X  10’  A/ cm’) 
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cycle  which  is  shown  in  Fig.  U.12.  The  figures  show  the  occurrence 
of  spatial  lag  similar  to  the  lag  seen  in  Fig.  4.11. 

Figures  U.17  and  4.18  show  the  time  variation  over  one  RF 
cycle  of  electric  field,  average  energy,  and  average  velocity  at 
fixed  points  inside  the  device.  The  locations  of  the  points  are 
shown  by  their  distances,  given  in  the  figures,  from  the  left- 
hand  contact.  One  point  is  located  near  each  of  the  two  contacts, 
and  a  third  near  the  metallurgical  junction.  The  energy  profiles 
show  temporal  lag  of  approximately  20  degrees  between  the  peak 
in  field  and  the  peak  in  majority  carrier  energy  at  the  points  near 
the  device  boxjndaries.  Wear  the  junction,  the  lag  is  somewhat  less. 
The  velocity  profiles  from  near  the  boundaries  show  undershoot  in 
majority  carrier  velocity  between  approximately  200-  and  2T0-degree 
phase,  and  overshoot  betweeen  approximately  300-  and  30-degree  phase. 
Similar  effects  occur  near  the  center  of  the  device,  but  to  a  lesser 
degree. 

The  phases  of  the  velocity  overshoots  and  undershoots  in 
Fig.  4.l8  correspond  to  the  phases  at  which  the  induced  current 
waveforms  in  Fig.  4.12  differ  in  shape.  It  appears  that  the  shape 
difference  in  the  induced  currents  is  due  to  undershoots  and  over¬ 
shoots  in  velocity,  while  the  lag  between  the  injected  and  the 
induced  current  waveforms  is  due  to  injection  delay  caused  by 
energy  lag.  The  shape  difference  is  roughly  symmetrical  about 
a  270-degree  phase,  so  it  probably  does  not  affect  negative  conduc¬ 
tance  strongly.  Evidently  the  increased  negative  conductance  seen 
in  energy  and  momentum  conserving  results  is  primarily  due  to  in¬ 
jection  delay  caused  by  energy  lag. 
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10’  A/cm’  AND  f  =  200  GHz) 


.  3  Effects  of  Boundary  Conditions 

Energy  and  momentum  consex^ing  simulation  involves  the  setting 
of  more  boundary  conditions  than  does  conventional  simulation,  so  it 
is  difficult  to  be  certain  when  conditions  used  intiie  two  types  of 
simulation  can  be  considered  equivalent.  The  boundary  conditions 
used  in  obtaining  the  energy  and  momentim  conserving  simulation 
results  presented  thus  far  were  those  described  in  Chapter  III.  Since 
the  inflow  boundary  conditions  may  influence  spatial  and  temporal 
lag,  thereby  affecting  device  admittance  results,  it  is  important 
that  their  influence  be  assessed. 

One  way  of  setting  realistic  boundai^  conditions  of  the  active 
region  of  a  device  is  to  add  a  highly  doped  contact  region  to  each 
end.  This  allows  the  energy  and  velocity  of  inflowing  carriers  to 
adjust  to  the  field  strength  in  the  contacts  before  the  carriers  enter 
the  active  region.  This  procedure  is  expensive  in  terms  of  computer 
time.  In  order  that  the  contact  regions  be  described  realistically, 
the  space  step  must  be  made  smaller  than  the  Debye  length  in  the 
contacts.  Numerical  stability  requires  that  At  be  reduced  in  pro¬ 
portion  to  Ax,  so  that  the  cost  of  simulation  goes  roughly  as  the 
square  of  the  number  of  space  steps. 

Some  simulation  of  the  0.3-um  device  with  contacts  added 
has  been  performed.  The  doping  profile  with  contacts  is  shown  in 
Fig.  U.19.  The  Debye  length  was  calcxilated  using^** 


fekgT 


(»*.l) 


where  N  is  the  doping  density  in  the  contacts.  At  a  temperature  of 
c 
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500°K,  the  Debye  length  in  the  contact  regions  shown  in  Fig.  it. 19 
is  approximately  37  X.  The  space  steps  and  time  steps  for  simulation 
of  the  device  with  contacts  were  chosen  to  be  26  %  and  O.OOU  ps, 
respectively.  This  time-step  length  is  short  enough  to  require  use 
of  a  Lax-Wendroff  type  finite-difference  scheme,  for  the  reasons 
discussed  in  Chapter  III.  Overall,  adding  the  contacts  shown  to  the 
0.3-um  device  increases  the  cost  of  simulation  by  more  than  a  factor 
of  ten. 

Figures  U.20  through  i*.23  show  profiles  of  electric  field, 
carrier  concentration,  average  energy,  and  average  velocity  at  four 
points  in  the  RF  cycle  when  the  0.3-pm  device  is  operated  at  300  GHz 
with  an  RF  amplitude  of  6  V.  Figures  U.2U  throu^  U.2T  show  similar 
results  with  no  contacts  present.  It  is  evident  that  the  behavior  of 
the  solution  of  the  transport  equations  across  the  active  region  is 
nearly  identical  in  the  two  cases.  It  may  be  concluded  that  it  makes 
little  difference  whether  specific  allowance  for  contact  regions  is 
made  or  whether  the  simple  boundary  conditions  described  in  Chapter 
III  are  applied. 

U.ii  Sources  of  Fhergy  Lag 

The  two  kinds  of  simulation  results  shown  in  Figs.  4.8  thro\igh 
4.10  diverge  more  rapidly  with  decreasing  device  length  than  with 
increasing  frequency.  This  can  be  seen  from  the  fact  that  100-GHz 
results,  which  have  been  obtained  for  each  of  the  three  devices, 
diverge  with  decreasing  length,  while,  for  a  given  device,  the 
difference  between  the  two  types  of  results  changes  little  with 
frequency.  The  way  in  which  the  results  behave  implies  that  the 
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FIG.  4.20  ELECTRIC  FIELD  (SOLID  LINES);  ELECTRON  DENSITY,  ELECTRON  ENERGY  AND  ELECTRON  VELOCITY 
(MINUS  SIGNS);  HOLE  DENSITY,  HOLE  ENERGY  AND  HOLE  VELOCITY  (PLUS  SIGNS)  VS.  DISTANCE 
IN  THE  0.3-Rin  DOUBLE-DRIFT  DEVICE  WITH  CONTACTS  AT  ZERO  DEGREES  PHASE  IN  THE  RF  CYCLE. 
{Vop  =  6  V,  J  =  1.5  X  10’  A/cjd"  and  f  =  300  GHz) 
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IN  THE  0.3-nni  DOUBLE-DRIFT  DEVICE  WITHOUT  CONTACTS  AT  ZERO  DEGREES  PHASE  IN  THE  RI  CYCLE. 
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FIG.  k.2J  ELECTRIC  FIELD  (SOLID  LINES);  ELECTRON  DENSITY^  ELECTRON  ENERGY  AND  ELECTRON  VELOCITY 
(MINUS  SIGNS);  HOLE  DENSITY,  HOLE  ENERGY  AND  HOLE  VELOCITY  (PLUS  SIGNS)  VS.  DISTANCE 
IN  THE  0.3-Mm  DOUBLE-DRIFT  DEVICE  WITHOUT  CONTACTS  AT  270  DEGREES  PHASE  IN  THE  RF  CYCLE. 


energy  lag  which  gives  rise  to  increased  negative  conductance  becomes 
more  pronoi.  ^ced  as  device  length  decreases,  but  changes  comparatively 
little  with  frequency.  Factors  which  might  contribute  to  the  length 
dependence  of  energy  lag  include  the  inflow  boundary  conditions  on 
energy,  carrier  cooling  due  to  Impact  ionization  by  opposing  carriers, 
and  spatial  variation  of  the  electric  field.  Each  of  these  will  now 
be  examined. 

Effects  of  the  inflow  boundary  conditions  on  spatial  lag  have 
been  tested  by  incorporating  "hot"  boundaries,  in  which  the  inflowing 
carriers  are  assigned  three  times  the  thermal  energy  associated  with 
the  lattice  temperature,  in  the  simulation  program.  This  might  be 
expected  to  reduce  the  amount  by  which  carriers  are  out  of  equili¬ 
brium  with  the  field  after  crossing  the  field  step  at  the  inflow 
bovindary.  Figures  U.28  and  U.29  show  resulting  profiles  of  electric 
field,  average  energy,  and  average  velocity  which  correspond  to 
those  shown  in  Figs.  U.2U  through  U.27.  The  similarity  between  the 
corresponding  profiles  shows  that  the  degree  of  energy  lag  has 
little  to  do  with  the  inflow  conditions  on  energy. 

Energy  lag  might  also  be  affected  by  the  cooling  of  each 
carrier  distribution  by  impact  ionizations  caused  by  the  opposing 
carrier  type.  This  possibility  will  be  explained  in  terms  of 
electrons.  The  number  of  electrons  entering  a  device  at  the  left- 
hand  boundary  is  small.  Just  upstream  of  the  botmdary,  impact 
ionizations  initiated  by  holes  may  produce  a  number  of  electrons 
comparable  to  or  greater  than  the  number  which  actually  cross  the 
boundary.  Consequently,  the  average  energy  of  electrons  near  the 
boimdary  may  be  determined  for  the  most  part  by  the  low  average 
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FIG.  4.28  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  ENERGY 

(MINUS  SIGNS)  AND  HOLE  ENERGY  (PLUS  SIGNS)  UNDER 
"HOT"  BOUT'DARY  CONDITIONS  VS.  DISTANCE  IN  THE 
0.3-mii  DOUBLE-DRIFT  DEVICE  AT  VARIOUS  POINTS  IN 
THE  RF  CYCLE.  (Vj^^  =  6  V,  =  1.5  x  10’  A/cm’ 
AND  f  =  300  GHe) 
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energy  of  those  electrons  which  are  produced  by  hole  ionizations. 

This  would  depress  the  average  energy  of  the  electron  distribution 
and  contribute  to  spatial  lag.  The  hole  distribution  near  the  right- 
hand  boundary  might  be  affected  in  a  similar  way. 

In  order  to  eliminate  this  "opposite-carrier  cooling"  effect, 
it  can  be  assumed  that  carriers  created  by  opposite-carrier  ioniza¬ 
tions  have  exactly  the  same  average  energy  as  those  already  present. 
Then  Eq.  2.32  would  become 

=  (w  -  w  )B  J  =  0  .  (U.2) 

Figures  i*.30  and  i*.31  show  electric  field,  average  energy,  and  average 
velocity  profiles  calculated  under  the  same  conditions  as  were  those 
shown  in  Fig.  k.2k  through  *+.27,  except  that  opposite-carrier  cooling 
has  been  eliminated  from  the  simulation  program.  The  change  clearly 
results  in  less  spatial  lag. 

These  results  suggest  a  reason  why,  as  seen  previously,  simula¬ 
tion  results  are  relatively  insensitive  to  changes  in  the  inflow 
boundary  conditions.  The  total  population  of  each  carrier  near  its 
inflow  boundary  conditions  is  dominated  by  carriers  produced  by 
opposite-carrier  ionizations.  Average  energy  and  velocity  near 
inflow  boundaries  are  determined  mainly  by  the  energy  and  velocity 
of  carriers  produced  by  impact  ionizations,  and  not  by  the  proper¬ 
ties  of  carriers  which  cross  the  inflow  boundary.  Thus  opposite- 
carrier  cooling  has  the  effect  of  decoupling  the  interior  of  the 
device  from  its  inflow  boundaries. 

Figures  U.32  and  U.33  show  admittance  results  from  the  1- 
and  0.3-pm  devices  with  and  without  opposite-carrier  cooling. 
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FIG.  Jt.30  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  ENERGY 

(MINUS  SIGNS)  AND  HOLE  ENERGY  (PLUS  SIGNS)  WITHOUT 
OPPOSITE-CARRIER  COOLING  VS.  DISTANCE  IN  THE  0.3-ym 
DOUBLE-DRIFT  DEVICE  AT  VARIOUS  POINTS  IN  THE  RF 
CYCLE.  (Vpp  =  6  V,  =  1.5  X  10*  A/cm’  AND 
f  =  300  GHz) 
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FIG.  4.31  ELECTRIC  FIELD  (SOLID  LINES),  ELECTRON  VELOCITY 

(MINUS  SIGNS)  AND  HOLE  VELOCITY  (PLUS  SIGNS)  WITHOUT 

OPPOSITE-CARRIER  COOLING  VS.  DISTANCE  IN  THE  0.3-vni 

DOUBLE-DRIFT  DEVICE  AT  VARIOUS  POINTS  IN  THE  RF 

CYCLE.  (V„  =  6  V,  =  1.5  X  10’  A/cm’  AND 
nr  uC 

f  =  300  GHz) 
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FIG.  4. 32  G-B  CHARACTERISTICS  FOR  THE  l-Miti,  DOUBLE-DRIFT 

DEVICE  WITH  (SOLID  LINE)  AND  WITHOUT  (DASHED  LINE) 
OPPOSITE-CARRIER  COOLING.  (V_^  =  10  V  AND 
J  .  =  6  X  10’  A/cin*) 


The  removal  of  cooling  from  the  simulation  can  be  seen  to  result 


in  increased  negative  conductance  In  both  devices,  with  the  increase 
being  much  more  pronounced  in  the  case  of  the  shorter  device. 

These  results  are  in  apparent  contradiction  to  the  idea  that  energy 
lag  is  what  causes  increased  negative  conductance  in  energy  and 
momentum  conserving  simulation  results.  Opposite-carrier  cooling 
increases  the  amoiont  cf  spatial  lag  associated  with  the  inflow  bound¬ 
aries,  but  it  reduces  negative  conductance.  This  can  be  seen  from 
comparison  of  Figs.  U.8  and  U.IO  with  Figs.  U.32  and  4.33.  Spatial 
lag  associated  with  opposite-carrier  cooling  clearly  does  not  cause 
the  observed  divergence  between  simulation  results.  Instead,  it 
tends  to  counteract  the  effects  of  whatever  does  give  rise  to  the 
divergence. 

The  mechanism  which  does  cause  the  two  types  of  results  to 
diverge  is  probably  associated  with  the  spatial  gradient  of  electric 
field  strength  inside  the  double-drift  devices.  In  contrast  to  the 
flat-field  structures  considered  in  Section  U.l,  the  double-drift 
structures  described  in  Section  4.2  are  doped  heavily  enough  to 
give  considerable  slope  to  the  field  strength.  The  slope  in  field 
gives  rise  to  a  kind  of  distributed  energy  lag  whose  degree  increases 
with  the  steepness  of  the  slope,  hence  with  increasing  doping  con¬ 
centration.  This  will  now  be  described  in  more  detail. 

The  total  time  rate  of  change  in  field  strength  seen  by  a 
moving  carrier  consists  not  only  of  the  time  rate  of  change  of 
field  at  the  current  position  of  the  carrier,  but  also  of  the  change 
seen  by  the  carrier  as  it  moves  through  spatial  variations  in  tile 
field.  The  total  rate  of  change  is  given  by 
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M  =  M  +  M 

dt  3t  ^  3x 


(U.3) 


During  much  of  the  buildup  of  the  injected  pulse,  the  carrier  current 
densities  are  small,  and  the  time  partial  of  the  electric  field  is 
approximately  equal  to  the  time  partial  of  the  terminal  voltage  divided 
by  the  device  length.  When  the  vbltage  is  sinusoidal,  the  maximum  of 
the  time  partial  of  the  field  at  angular  frequency  lo  is  given  by 


3t  L 


(U.l*) 


and  occurs  at  0-degree  phase.  When  the  carrier  density  is  small,  the 
space  partial  of  the  field  is  related  to  the  doping  density  by 
Poisson's  equation: 


3E 

3x 


(h.5) 


For  electrons  in  the  p  layer  of  the  1-ym  double-drift  device 
at  a  frequency  of  100  GHz  and  an  RF  amplitude  of  10  V,  Eq.  4.3  becomes 


—  =  6.3  X  10^®  V/cm*s  +  1.2  x  10^^  V/cm»s  ,  (4.6) 

dt 

where  it  is  assumed  that  u  is  8.5  x  10*  cm/s.  Similarly,  in  the 
0.3-wm  device  at  an  amplitude  of  6  V,  the  equation  becomes 


=  1.3  X  10^^  V/cm*s  +  2.9  x  10^’  V/cm*8  .  (4.7) 

dt 

Equations  4.6  and  4.7  show  that  the  total  time  rate  of  change 
of  field  seen  by  a  moving  carrier  at  a  given  frequency  is  much  larger 
in  the  0.3-wm  device  than  in  the  1-pm  device.  Hie  equations  also 
show  that,  even  at  the  moment  when  the  time  partial  of  the  field  is 
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at  its  maximum,  the  spatial  gradient  of  field  contributes  the 
majority  of  the  total  rate  of  change  seen  by  a  moving  carrier. 

The  greater  the  total  rate  of  change,  the  greater  the  lag  between 
carrier  energy  and  field,  so  that  at  the  center  of  the  0.3-pm 
device  the  ionization  rates  peak  later  in  the  RF  cycle  than  in  the 
1-ym  device,  giving  rise  to  more  injection  delay  and  a  larger  i.icrease 
in  negative  conductance  relative  to  the  conventional  result. 

k.'y  Limitations  on  IMPATT  Performance 

The  simulation  results  presented  in  this  chapter  indicate  that 
Si  diodes  will  support  operation  of  the  IMPATT  mode  at  frequencies  up 
to  at  least  300  GHz.  It  will  be  shown  in  this  section  that  the 
material  properties  do  eventually  impose  a  fundamental  upper  frequency 
limit  on  IMPATT  operation.  The  frequency  at  which  this  limit  lies 
can  be  estimated  for  any  material  by  use  of  the  energy  balance  rela¬ 
tion  and  is  in  the  submillimeter-wave  region  for  Si  devices.  This 
section  also  discusses  why  the  simulation  predicts  millimeter-wave 
device  efficiencies  which  are  in  excess  of  those  which  are  obtained 
in  experiment.  It  is  shown  that  this  can  be  fully  explained  by  the 
presence  of  parasitic  series  resistance  external  to  the  active 
IMPATT  layer. 

The  existence  of  an  upper  frequency  limit  for  IMPATT  operation 
is  a  result  of  the  way  in  which  carrier  energy  responds  to  the  time 
variation  of  the  electric  field.  The  simple  energy  balance  rela¬ 
tion,  which  accounts  for  energy  gain  from  the  field  and  loss  to 
lattice  collisions,  provides  an  approximate  description  of  the 
energy  response; 
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dv 

dt 


w  -  w 


ik.8) 


=  quE  -  — - — -  . 

w 

If  u  is  constant  and  the  time-varying  component  of  the  field  is  given 
by  the  real  part  of  the  solution  for  the  time-varying  component 

of  energy  is 

(U.9) 

The  energy  response  given  by  Eq.  k.9  follows  the  familiar 

single-pole  transfer  function.  A  Bode  plot^^  for  the  normalized 

response  is  shown  in  Fig.  U.3l*.  The  plot  shows  that  the  amplitude  of 

the  response  is  down  by  a  factor  of  two  at  a  frequency  of  1/2itt  . 

w 

In  view  of  the  rapid  variation  of  the  ionization  rates  with  energy 
which  is  shown  in  Fig.  2.3,  it  appears  that  IMPATT  mode  operation  will 
be  seriously  degraded  at  this  frequency  because  of  reduced  particle 
current  modulation.  The  response  begins  to  roll  of f  an  octave  lower, 
implying  that  material  properties  (specifically  the  energy  relaxa¬ 
tion  time)  will  impose  an  upper  frequency  limit  of  approximately 

1/4itt  for  substantially  undegraded  IMPATT  operation, 
w 

Figure  4.3^  also  shows  the  increasing  phase  lag  between 
energy  and  field  as  frequency  increases.  Lippens  and  Constant**® 
have  used  the  energy  balance  relation  to  support  previously  pub¬ 
lished  findings  from  this  work  which  pointed  to  energy  lag  as  the 
reason  for  predictions  of  increased  negative  conducteuice  in  energy 
and  momentum  conserving  simulation  results.  Lippens  and  Constant 
suggested  that  the  findings  could  be  explained  in  terms  of  the  phase 
lag  shown  in  the  figure.  While  their  conclusion  that  the  phase  lag 
in  itself  tends  to  make  IMPATT  operation  more  efficient  is  valid. 


w  =  Re 


w  „  amt 

- qE  ue 

1  +  1U)T  o 


-124- 


they  overlooked  the  roll  off  in  the  energy  response  which  is 
predicted  by  the  energy  balance  relation.  It  should  be  noted 
that  the  simple  balance  relation  fails  to  allow  for  spatial  inhomo¬ 
geneities,  In  fact,  as  previously  noted,  spatial  field  gradients 
are  apparently  the  principal  reason  for  differences  between 
conventional  and  energy  and  momentum  conserving  simulation  results 
for  submicron  IMPATT  devices. 

While  amalysis  involving  the  energy  balance  relation  is 

approximate,  a  more  detailed  analysis  would  not  invalidate  the 

conclusion  that  the  energy  response  will  roll  off  at  high  frequencies. 

Equation  U.9  predicts  that  the  frequency  at  which  roll  off  begins 

is  approximately  l/4irT  .  The  energy  relaxation  time  for  electrons 

w 

in  Si  was  estimated  in  Chapter  II  to  be  approximately  0.15  ps 
over  a  broad  range  of  energy,  implying  that  IMPATT  mode  operation 
in  Si  will  deteriorate  substantially  at  frequencies  above  500  GHz. 
This  roll  off  in  the  energy  response  may  be  the  reason  why  the  two 
admittance  curves  in  Fig.  4.10  cross  Just  below  500  GHz. 

The  best  reported  experimental  efficiencies  for  millimeter- 
wave  Si  IMPATTs  are  1  percent  or  less  at  frequencies  above  150 
GHz.’*  This  is  much  lower  than  the  efficiencies  routinely  obtained 
from  microwave  devices.  The  results  of  this  study  suggest  that  this 
is  not  due  to  an  intrinsic  failure  of  the  IMPATT  mode  in  Si  at 
millimeter-wave  frequencies.  An  alternative  explanation  is  that 
parasitic  losses  are  responsible  for  low  millimeter-wave  efficien¬ 
cies  and  constitute  the  major  limitation  on  the  experimental  per¬ 
formance  of  millimeter-wave  IMPATT  diodes.  Sources  of  parasitic 
lose  include  undepleted  substrate  layers  in  diode  structures, 


-126- 


ohmic  contact  resistance,  series  resistance  of  bond  wires,  and  other 
circuit  losses.  Loss  mechanisms  which  are  unimportant  at  lower 
frequencies  may  become  important  in  the  millimeter-wave  range  be¬ 
cause  of  skin  effects,  surface  scattering  effects,  and  reductions 
in  carrier  mobility  due  to  velocity  relaxation. 

A  detailed  analysis  of  parasitic  loss  is  outside  the  scope 
of  this  study,  but  the  effect  of  loss  cam  be  estimated  by  placing 
an  equivalent  parasitic  resistance  in  series  with  the  "intrinsic" 
diode.  Figure  4.3^  shows  both  "intrinsic"  negative  resistance  and 
efficiency  in  the  presence  of  various  values  of  series  resistamce 
as  functions  of  RF  aimplitude  for  the  0.3  pm  device  operating  at 
300  GHz.  A  device  diameter  of  0.5  mil  is  assumed.  "Intrinsic" 
efficiency  as  predicted  by  conventional  simulation  is  also  shown. 

The  peak  efficiency  predicted  by  energy  and  momentum  conserving 
simulation  is  much  greater  than  that  predicted  by  conventional 
simulation,  and  both  efficiencies  are  large  enough  to  indicate 
little  deterioration  in  the  operation  of  the  IMPATT  mode  in  Si 
at  frequencies  up  to  300  GHz.  The  presence  of  parasitic  resistance 
can  be  seen  to  reduce  efficiency  drastically. 

4.6  Summary  and  Conclusions 

Simulation  results  show  that  the  additional  physical  effects 
allowed  for  in  the  energy  and  momentum  conserving  transport  model 
cause  the  average  energy  of  carriers  to  lag  behind  the  local  electric 
field.  Energy  lag  can  occur  in  space,  downstream  of  abrupt  changes 
in  field,  or  it  can  occur  in  time,  in  the  presence  of  rapid  time 
variation  of  field.  Energy  lag  gives  rise  to  overshoots  and  under¬ 
shoots  in  carrier  velocity  as  compared  to  conventional  field-dependent 
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FIG.  4.35  "INTRINSIC"  NEGATIVE  RESISTANCE  (DASHED  LINE) 

AND  EFFICIENCY  IN  THE  PRESENCE  OF  VARIOUS  VALUES 
OF  SERIES  RESISTANCE  (SOLID  LINES)  VS.  RF  AMPLITUDE 
FOR  THE  0.3-vm  DOUBLE-DRIFT  DEVICE,  (f  =  300  GHz 


AND  J 


1.5  X  10’  A/cjb*) 


drift  velocity.  Local  departures  from  the  conventional  velocity 
can  be  large,  but  tend  to  affect  minority  carriers  rather  than 
majority  carriers  in  IMPATTs,  so  their  effect  on  device  terminal 
admittance  is  generally  small. 

As  the  double-drift  device  length  becomes  less  than  O.5  Mm, 
energy  and  momentum  conserving  simulation  predicts  substantially 
better  IMPATT  performance  than  does  conventional  simulation.  A 
70-percent  difference  in  predicted  optimum  efficiency  at  300  GHz 
has  been  observed.  Differences  in  predicted  performance  appear  to 
be  due  to  injection  delay  caused  by  the  lag  between  the  carrier 
energy  and  the  electric  field.  This  lag  appears  to  be  caused  by 
spatial  field  gradients. 

Simulation  results  are  relatively  insensitive  to  boundary 
conditions  on  carrier  energy,  and  the  boundary  conditions  described 
in  Chapter  II  have  been  found  to  give  results  consistent  with  those 
obtained  incorporating  realistic  contact  regions.  The  degree  of 
carrier  cooling  which  results  from  impact  ionization  by  opposing 
carriers  is  of  more  importance  to  device  behavior  than  boundary 
conditions.  This  cooling  increases  the  amount  of  spatial  energy 
lag  and  reduces  negative  conductance. 

Simulation  results  predict  greater  millimeter-wave  device 
efficiencies  than  are  observed  in  experiment.  This  can  be  accounted 
for  by  the  presence  of  modest  amounts  of  parasitic  series  resistance. 
It  Is  apparently  such  resistance,  rather  than  a  deterioration  of 
IMPATT  mode  operation  at  high  millimeter-wave  frequencies,  that 
limits  the  performance  of  experimental  devices.  In  the  absence 
of  parasitic  loss,  the  upper  frequency  of  operation  of  the  IMPATT 
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mode  would  be  subject  to  a  fundamental  limit  which  is  set  by 


material  properties.  This  limit  is  estimated  to  be  approximately 
500  GHz  for  Si  devices. 
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CHAPTER  V.  SUMMARY,  CONCLUSIONS, 

AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

^.1  Sununary  and  Conclusions 

The  purpose  of  this  study  was  to  develop  and  apply  a  new 
class  of  semiconductor  device  simulation  which  is  more  general  than 
existing  drift-diffusion  based  simulations,  and  to  use  the  Simula- 
tion  to  investigate  the  "intrinsic"  properties  of  millimeter- 
'  wave  Si  IMPATTs.  A  self-consistent  bipolar  energy  and  momentum 

conserving  simulation  has  been  developed  for  this  purpose.  The 
work  that  was  performed  falls  into  three  categories:  model  defini¬ 
tion,  model  implementation,  and  model  utilization.  The  overall 
effort  is  believed  to  be  the  first  reported  application  of  a  self- 
consistent  bipolar  energy  and  momentum  conserving  transport  model 
to  semiconductor  device  simulation. 

Work  performed  in  the  category  of  model  definition  was 
described  in  Chapter  II  of  this  dissertation.  The  charge-transport 
model  developed  consists  of  transport  equations  for  electron  and 
•  hole  concentration,  average  velocity,  and  average  energy  as 

functions  of  time  and  space.  It  provides  a  second-order  descrip- 
tion  of  the  carrier-velocity  distributions,  and  allows  the  dis¬ 
tributions  to  evolve  in  space  and  time  in  a  self-consistent  manner. 
This  is  in  contrast  to  the  conventional  model  based  on  drift  and 
diffusion,  which  implicitly  assumes  that  the  distributions  are  always 
in  equilibrium  with  the  local  electric  field. 
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Ihe  transport  equations  themselves  are  the  zeroth,  first,  and 


second  velocity  moments  of  the  phase-space  transport  equation.  Except 
for  the  terms  describing  the  effects  of  collisions,  the  equations  were 
obtained  by  application  of  the  method  of  moments.  The  collision  terms 
were  not  obtained  in  this  way  because  the  form  of  the  phase-space 
collision  term  for  carriers  in  Si  is  not  sufficiently  well  known 
to  permit  the  taking  of  velocity  moments.  Instead,  forms  for  the 
collision  terms  were  chosen  in  such  a  way  as  to  allow  for  their 
dependence  on  the  concentration,  average  energy,  and  average 
velocity.  The  collision  parameters  used  in  the  model  are  energy- 
dependent  relaxation  times  for  energy  and  velocity,  and  energy- 
dependent,  per-unit-time  impact  ionization  rates.  The  latter  are 
more  ftindamental  than  the  conventional  per-vinit-di stance  ionization 
rates  because  they  do  not  assiome  any  correlation  between  the  average 
carrier  energy,  which  determines  the  number  of  ionizations  taking 
place  per  unit  time,  and  the  average  velocity,  which  determines  the 
distance  traveled  in  the  average  time  between  ionizations. 

The  parameters  were  given  numerical  values  by  simplifying 
the  transport  equations  to  describe  the  situation  of  a  spatially 
uniform,  dc  electric  field,  and  requiring  that  results  from  the 
simplified  equations  be  consistent  with  experimentally  known  values 
of  electron  and  hole  drift  velocities  and  ionization  rates.  Mapping 
the  parameters  onto  the  carrier  energies  was  accomplished  through 
use  of  a  theoretically  determined  equilibrium  relationship  between 
the  carrier  energies  and  the  electric  field.  The  hole  and  electron 
energy  relaxation  times  determined  in  this  way  were  approximately 
0.075  and  0.15  ps,  respectively,  over  much  of  the  carrier  energy 
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range  considered.  The  velocity  relaxation  times  were  fovmd  to 
rainge  from  twice  the  energy  relaxation  times  at  low  energies  to 
less  than  0.02  ps  at  high  energies. 

A  correspondence  between  the  energy  and  momentum  conseirving 
transport  model  and  the  conventional  drift-diffusion  model  was 
established.  This  showed  that,  under  the  conditions  of  slow  space 
and  time  variation  of  the  electric  field,  the  conventional  diffu¬ 
sion  coefficient  can  be  written  in  terms  of  the  velocity  relaxation 
time  and  the  temperature  of  the  carrier  distribution.  The  result¬ 
ing  estimates  for  the  electron  and  hole  diffusion  coefficients  as 
functions  of  the  electric  field  drop  off  considerably  from  their 
low-field  values  as  field  strength  increases.  The  energy  and 
momentm  conserving  model  was  also  compared  to  other  nonconventional 
models  which  have  been  applied  to  modeling  of  IMPATTs.  The  latter 
were  shown  to  be  less  complete  and  less  self-consittent  than  the 
energy  and  momentum  conserving  model. 

Computer  implementation  of  the  transport  model  using  finite- 
difference  techniques  was  described  in  Chapter  in,  A  set  of 
normalizations  was  developed  which  (for  constant  time  and  space 
step)  eliminates  the  time  step,  the  space  step,  and  several  phy¬ 
sical  constants  from  the  finite-difference  equations.  Stability 
analysis  of  various  finite-difference  schemes  showed  that  many 
explicit  forms  of  the  energy  and  momentum  conserving  transport 
equations  are  unstable  when  the  time  step  is  short  in  comparison 
to  the  relaxation  times,  and  other  stable  forms  may  exhibit 
undesirably  large  amounts  of  numerical  diffusion.  A  scheme  of 
the  L8uc-Wendroff  type  was  found  to  be  capable  of  giving  stability 
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and  minimal  numerical  diffusion  for  eO.!  time-step  lengths.  It 
was  further  shown  that,  when  the  time  step  is  sufficiently  long, 
certain  more  efficient  schemes  which  are  otherwise  unstable  can  be 
made  stable  by  the  use  of  an  advanced-time  form  of  the  energy 
and  velocity  relaxation  terms. 

Chapter  IV  presented  results  of  computer  simulation  of 
millimeter-wave  Si  avalanche  diode  structures.  Simulation  of  low- 
doped  diodes  with  spatially  \miform  electric  field  shows  that 
rapid  field  changes  in  space  (such  as  at  the  boundaries  of  deple¬ 
tion  regions)  and  time(  as  in  the  presence  of  large-signal  RF 
terminal  voltages)  can  give  rise  to  nonequilibrium  between  the 
carrier  velocity  distribution  and  the  electric  field.  Such  non¬ 
equilibrium  manifests  itself  in  the  form  of  lag  between  the  carrier 
energy  and  local  electric  field,  and  gives  rise  to  velocity  over¬ 
shoot  or  undershoot.  The  length  of  the  distance  over  which  the 
minority  carrier  distribution  adjusts  itself  to  a  field  step  in 
space  was  shown  to  be  strongly  affected  by  the  initial  temperature 
assigned  to  carriers  generated  by  impact  ionization  by  carriers  of 
the  "opposing"  type. 

Qiergy  and  momentum  conserving  simulation  results  were  found 
to  diverge  significantly  from  conventional  results  as  device  length 
becomes  less  than  1  pm.  Energy  and  momentum  conserving  results 
predict  greater  negative  conductance  than  conventional  results, 
apparently  because  of  Increased  Injection  delay  caused  by  lag 
between  the  carrier  energies  and  the  electric  field.  The  amount 
by  which  the  two  types  of  results  diverge  depends  more  strongly 
on  device  length  than  upon  operating  frequency.  This  Is  apparehtly 


due  to  the  increased  field  gradients  which  occur  as  device  length 
decreases.  Differences  between  predicted  performances  at  optimum 
RF  amplitude  are  particularly  large. 

Simulation  results  showed  no  degradation  in  operation  of  the  IMPATT 
mode  in  Si  at  frequencies  up  to  300  GHz.  Calculations  based  on  the 
energy  balance  relation  predicted  that  the  IMPATT  mode  in  Si  will 
begin  to  degrade  as  frequency  exceeds  500  GHz  because  the  carrier 
energies  will  cease  to  respond  to  time  variation  of  the  electric 
field.  The  millimeter-wave  efficiencies  predicted  by  simulation 
results  are  well  in  excess  of  those  which  have  been  obtained  in 
experiment.  This  strongly  suggests  that  parasitic  loss,  rather  than 
failure  of  the  intrinsic  IMPATT  mode,  is  the  factor  which  presently 
limits  the  performance  of  the  current  state-of-the-art  devices. 

5.2  Suggestions  for  Further  Research 

The  present  investigation  has  laid  a  foundation  for  a  variety 
of  further  work.  Suggestions  for  further  research  can  be  considered 
in  three  categories:  transport  modeling,  numerical  methods,  and 
device  simulation. 

The  most  obvious  extension  in  the  category  of  transport  model¬ 
ing  is  to  incorporate  nonequivalent  conduction  band  valleys.  This 
would  extend  the  applicability  of  the  simulation  to  III-V  comi>o\inds 
such  as  GaAs  and  InP,  which  are  widely  used  to  fabricate  high- 
performance  semiconductor  devices.  It  is  known  that  certain  relaxa¬ 
tion  times  in  GaAs  cuid  InP  are  much  longer  than  those  encountered 
in  Si,’’  and  so  "transient"  effects  will  probably  be  of  Importance 
at  lower  frequencies  and  in  larger  devices  than  is  the  case  f6r  Si. 
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TSie  forms  of  the  equations  for  the  multi-valley  case  have  been 
presented  by  Blotekjaer. Hovever,  compared  to  the  single¬ 
valley  situation,  additional  relaxation  times  are  needed  to  char¬ 
acterize  intervalley  transfer.  These  cannot  be  determined  by  the 
method  used  in  Chapter  II  (i.e.,  requiring  consistency  vith 
experimental  measurements),  and,  as  discussed  below,  will  have  to 
be  obtained  in  some  other  way.  Work  on  extending  the  energy  and 
momentum  conserving  model  to  transport  in  GaAs  is  currently  underway 
within  the  Electron  Physics  Laboratory. 

Under  low-field  conditions  calculation  of  the  required 
additional  relaxation  times  is  reasonably  straightforward,  using 
statistical  Monte  Carlo  simulation  or  the  Rees  iterative  technique  .** "* 
Existing  theoretical  calculations  seem  to  be  inconsistent  with 
experimental  observations  for  intermediate  «uid  high  fields — 
typically  above  20  kV/cm.  (One  example  is  the  calculation  by 
Jacoboni  et  al.^®  of  mean  carrier  energies  of  the  order  of  0.5  eV 
at  a  field  of  100  kV/cm.  This  implies  significant  impact  ionization 
at  this  field,  something  which  is  not  observed  in  practice.)  The 
discrepancy  between  theory  and  experiment  may  be  due  to  neglect  of 
the  intra-collisional  field  effect.”  A  major  priority  for  future 
transport  modeling  should  be  transport  characterization  including 
impact  ionization  and  the  intra-collisional  field  effect.  The 
deterministic  Rees  iterative  technique  is  probably  best  suited  for 
transport  characterization  and  parameter  evaluation,  since  it 
avoids  the  large  statistical  variance  associated  with  Monte  Carlo 
characterization  of  "rare”  events  such  as  impact  ionization,  and 


is  also  well  suited  to  describing  nonlinear  effects  such  as 
transport  in  degenerate  semiconductors  and  Auger  recombination. 

Additional  work  in  the  area  of  numerical  methods  shoiild 
include  the  extension  to  two  spatial  dimensions.  This  will  probably 
not  present  fundamental  problems;  the  stability  criteria  for  finite- 
difference  approximations  to  the  two-dimensional  transport  equations 
should  be  similar  to  those  established  in  this  work  for  the  one¬ 
dimensional  case,  and  techniques  for  the  rapid  solution  of  the  two- 
dimensional  Poisson's  equation  are  now  well  established. ® “  The 
extension  to  two  dimensions  would  permit  energy  and  momentum  con¬ 
serving  simulation  of  transistor  structures.  Another  extension  of 
the  simulation  would  be  incorporation  of  a  nonzero  heat  flow 
vector.  It  does  not  appear  that  this  would  be  difficult,  although 
stability  of  proposed  finite-difference  equations  incorporating  a 
description  of  heat  flow  should  be  examined  using  the  techniques 
applied  in  Chapter  III.  Finally,  renewed  attention  could  be  given 
to  the  prospects  of  obtaining  an  implicit  or  semi-implicit  numerical 
scheme  whose  stability  would  permit  use  of  longer  time  steps.  The 
advantage  of  such  a  scheme  is  not  certain,  however,  since  time-step 
size  is  limited  by  accxxracy  requirements  even  for  unconditionally 
stable  schemes. 

The  application  of  energy  and  momentum  conserving  simulation 
to  semiconductor  device  modeling  is  in  its  infancy.  Advances  in 
technology  have  only  recently  made  possible  the  fabrication  of 
devices  so  small  that  traditional  drift-diffusion  simulation  pro¬ 
vides  a  generally  Inappropriate  description.  The  overshoot 
effects  in  small  devices  can  be  described  using  parameterized 


distributions  as  in  the  present  work,  or  using  "exact" 
distributions.  The  former  approach  is  much  less  expensive. 

Only  when  a  case  can  be  made  that  the  "fine  structure"  of  the 
distribution  is  significant  (as  in  the  Jones-Rees  effect**  in  Gunn 
diodes)  does  it  seem  worthwhile  attempting  to  use  "exact"  distri¬ 
bution  models.  The  previously  mentioned  extensions  to  the  present 
work  would  yield  an  energy  and  momentum  conserving  simulation 
applicable  to  a  wide  range  of  semiconductor  devices,  and  could  provide 
answers  to  many  questions  of  current  interest  in  the  field  of  sub¬ 
micron  electronics.  Possible  one-dimensional  simulation  investi¬ 
gations  include  simulation  of  other  tr2msit-time  devices  (including 
TUNNETTs  and  transferred  electron  devices)  and  investigation  of 
the  large-signal  RF  properties  of  nn^  junctions.  A  two-dimensional 
simulation  could  be  used  to  model  short-channel  effects  in  MESFETs 
and  MOSFETs,  and  to  examine  the  feasibility  of  proposed  "novel" 
transistors  such  as  the  ballistic  transistor**  and  planar  doped 
barrier  transistor.*’  Overall,  energy  and  momentum  conserving 
simulation  seems  likely  to  become  the  standard  approach  of  sub¬ 
micron  semiconductor  device  modelers. 

The  most  significant  mesuis  of  comparison  between  conventional 
and  energy  and  momentum  conserving  IMPATT  simulation  results  may 
be  to  compare  on  the  basis  of  maximum  obtainable  efficiency.  A 
logical  extension  of  the  work  presented  in  this  study  would  be  to 
use  the  two  types  of  simulation  to  search  for  optimum  IMPATT 
structures  and  large-signal  operating  points  as  functions  of 
frequency  over  the  millimeter-wave  range. 


APPENDIX  A.  MATERIAL  PARAMETERS 


The  material  parameters  used  throughout  this  study  pertain 
to  Si  at  a  lattice  temperature  of  500°K.  The  diffusion  coefficients 
used  in  all  drift-diffusion  simulation  are  those  given  in  Fig.  2.4. 
Static  ionization  rates*'  are  given  in  Table  A.l.  Table  A. 2  gives 
static  drift  velocities,  and  Table  A. 3  lists  a  number  of  other 
parameters . 


Table  A.l 

Static  Ionization  Rates*' 
a(E)  -  A  exp  (-  b/E)  cm”^ 


Quant 5  ty 

Holes 

Electrons 

E  (kV/cm) 

A  (cm~^) 

2.0  X  10® 

2.6  X  10® 

0  <  E  <  240 

2.0  X  10® 

6.2  X  10* 

240  <  E  <  530 

5.6  X  10* 

5.0  X  10* 

E  >  530 

b  (kV/cm) 

2.17  X  10® 

1.69  X  10® 

0  <  E  <  24o 

2.17  X  10® 

1.31  X  10® 

240  <  E  <  530 

1.54  X  10® 

1.25  X  10® 

E  >  530 

Table  A. 2 

Static  Drift  Velocity® 

Quantity  Holes  Electrons 

Vgat  (cm/s)  1.02  x  lo’  8.5  x  10* 

W  (cm*/V*s)  250  550 
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Table  A. 3 


Other  Parameters 


Quantity 

Symbol 

Value 

Dielectric  constant 

e 

l.Oli  X  10"^^  F/cm 

Ionization  threshold 
energy 

E 

c 

2.0  eV 

Optical  phonon 
energy 

7[ui 

0.056  eV 

Effective  mass 

m 

**.55 

8.8U 

X  10”®^  kg  (holes) 

X  10~*^  kg  (electrons 

Mean  free  path 

X 

60  X  (holes) 

80  X  (electrons) 
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